143 Commits

Author SHA1 Message Date
Sven Czarnian
e19e50b027 update to the current interface version 2021-12-29 09:37:23 +01:00
Sven Czarnian
c5d5e332c8 fix compiler issues 2021-12-26 10:59:22 +01:00
Sven Czarnian
7411955ede fix a typo 2021-12-26 10:44:05 +01:00
Sven Czarnian
1873c68ab0 use the server's IAFs to check if we have a new inbound and over which IAF 2021-12-26 10:43:15 +01:00
Sven Czarnian
0a586f88c3 remove double lock 2021-12-26 09:29:50 +01:00
Sven Czarnian
06271fde5b code formatting 2021-12-26 09:29:37 +01:00
Sven Czarnian
4e635514f1 initialize the complete inbound after the creation 2021-12-23 17:50:35 +01:00
Sven Czarnian
d24c1ee1ce add a flag to indicate if the filter needs to be reset 2021-12-23 17:50:21 +01:00
Sven Czarnian
f54711f490 add a missing initialization 2021-12-23 17:50:08 +01:00
Sven Czarnian
031c9faa97 add a low pass filter to keep the TTL more stable 2021-12-23 17:49:56 +01:00
Sven Czarnian
d7c41d1941 give the controller the information he needs and that he can configure the tags 2021-12-23 09:58:31 +01:00
Sven Czarnian
ec9e9e285b predict the planned waypoint if a heading is given 2021-12-20 11:08:07 +01:00
Sven Czarnian
2795011021 check if the PTA exists 2021-12-19 13:06:19 +01:00
Sven Czarnian
c3a75bb900 add a menu to force set runways 2021-12-19 12:23:23 +01:00
Sven Czarnian
c14698bc39 fix the WTC message 2021-12-17 13:32:30 +01:00
Sven Czarnian
4f84e9a051 remove the debug message 2021-12-17 13:05:38 +01:00
Sven Czarnian
861232d09d update the IAF distance estimation to provide better results 2021-12-16 16:02:35 +01:00
Sven Czarnian
791c9c2605 do not filter based on distances 2021-12-16 16:02:08 +01:00
Sven Czarnian
4da604af99 add a function to project the IAF on the current path to have a virtual IAF 2021-12-16 16:01:49 +01:00
Sven Czarnian
ffb944fd6d add the remaining distance 2021-12-16 16:01:05 +01:00
Sven Czarnian
7238d6b92e use the correct return values 2021-12-16 16:00:44 +01:00
Sven Czarnian
1570b9c7f1 make the matcher public 2021-12-16 16:00:31 +01:00
Sven Czarnian
85600fcd57 add Eigen as a dependency 2021-12-16 16:00:12 +01:00
Sven Czarnian
27b422d806 add Eigen as a dependency 2021-12-16 15:59:56 +01:00
Sven Czarnian
a08ac9b244 fix compiler issues 2021-12-15 16:55:07 +01:00
Sven Czarnian
dbb3d6d557 fix the issues with the predicted path to planned path matching 2021-12-15 16:53:55 +01:00
Sven Czarnian
88c40575ff do not delete the airports 2021-12-15 13:25:04 +01:00
Sven Czarnian
de25092e3b fix sweatbox issues 2021-12-15 13:24:55 +01:00
Sven Czarnian
c94847b607 allow more tag entries to customize the views 2021-12-15 13:24:42 +01:00
Sven Czarnian
e108d3de57 fix some issues with the path prediction and add some addional information 2021-12-15 13:23:33 +01:00
Sven Czarnian
f053ead918 send the new sequence after every update 2021-12-09 12:29:38 +01:00
Sven Czarnian
ef91a3badd cleanup internal structures as soon as the controller is disconnected 2021-12-09 12:29:24 +01:00
Sven Czarnian
b88d2aabda set the airport 2021-12-09 10:27:40 +01:00
Sven Czarnian
076109be41 fix a return value bug 2021-12-09 10:27:30 +01:00
Sven Czarnian
0fc113cd2d use the standard color and a * to indicate if a plan is fixed or not 2021-12-02 19:07:03 +01:00
Sven Czarnian
c6159144bb indicate if the reload was successful 2021-12-02 19:06:23 +01:00
Sven Czarnian
6152c40dca fix compiler issues 2021-12-02 19:06:13 +01:00
Sven Czarnian
95493674e8 add the ignore for used parameters 2021-12-02 19:06:03 +01:00
Sven Czarnian
ae0442b560 add a menu 2021-12-02 19:05:34 +01:00
Sven Czarnian
41a9fd20d5 fix code analysis findings 2021-11-30 08:50:28 +01:00
Sven Czarnian
e4ea8f46c8 ignore the sweatbox and playbacks to use only real VATSIM data 2021-11-29 21:00:37 +01:00
Sven Czarnian
055cf900fb fix some changes and add a low pass filter to avoid fast ETA changes 2021-11-29 19:21:29 +01:00
Sven Czarnian
d89e8fef4c remove a debug message 2021-11-29 19:20:41 +01:00
Sven Czarnian
684ca40c1c fix a typo in the scaling 2021-11-25 22:51:25 +01:00
Sven Czarnian
eb521fca8b update the inbound queue lock 2021-11-25 22:28:32 +01:00
Sven Czarnian
1936b3649b use one second resolution 2021-11-25 22:27:34 +01:00
Sven Czarnian
8390cd5309 define the inbound and predict the flight path 2021-11-25 22:27:03 +01:00
Sven Czarnian
1594be8dc1 adapt to avoid blocking calls 2021-11-25 22:26:26 +01:00
Sven Czarnian
f2f68708ac register the function and the tag 2021-11-25 22:26:11 +01:00
Sven Czarnian
9146446924 change the API 2021-11-25 22:25:51 +01:00
Sven Czarnian
874cc7871d add getter for the arrival waypoints 2021-11-25 22:25:11 +01:00
Sven Czarnian
2240978b01 define the arrival waypoint 2021-11-22 16:18:47 +01:00
Sven Czarnian
155a035a9c define the wind data 2021-11-22 16:18:29 +01:00
Sven Czarnian
319e5b2952 build the geographic coordinates 2021-11-22 16:18:16 +01:00
Sven Czarnian
92bf63bf45 add the helper 2021-11-22 16:17:35 +01:00
Sven Czarnian
969f448bf4 define the geographic library 2021-11-22 16:17:20 +01:00
Sven Czarnian
9f80ce1631 add the helper for the UTC time and mach checks 2021-11-22 16:17:07 +01:00
Sven Czarnian
09e29afe7b add GeographicLib 2021-11-22 16:16:36 +01:00
Sven Czarnian
b1b6b5c790 do not send empty lists 2021-11-21 08:54:22 +01:00
Sven Czarnian
ce23f956ed code formatting 2021-11-21 08:54:11 +01:00
Sven Czarnian
b309e86ec6 change the configuration to handle the new ZMQ-REQ-REP-protocol and the RestAPI-Update protocol 2021-11-21 08:53:58 +01:00
Sven Czarnian
0acb45dd27 force a reload 2021-11-20 09:29:01 +01:00
Sven Czarnian
1bb9813a7a allow to add flights as forced for AMAN to include them, if the neighboring controller is not using AMAN 2021-11-20 09:28:36 +01:00
Sven Czarnian
2cc7c401e8 add the stands of GrPlugin and TST 2021-11-18 16:02:52 +01:00
Sven Czarnian
689c2326c2 use only one backend and adapt to the new communication structure
avoid empty messages
2021-11-18 16:02:41 +01:00
Sven Czarnian
a65aa445fa compare the backend's version with the plugin's version 2021-11-18 15:55:37 +01:00
Sven Czarnian
25c1519bb8 update to the current version 2021-11-18 15:53:57 +01:00
Sven Czarnian
a23de5ee67 add required libraries 2021-11-18 15:53:34 +01:00
Sven Czarnian
1a58301d4f send the reports every ten seconds 2021-11-13 23:03:34 +01:00
Sven Czarnian
d3551c9857 fix a compiler warning 2021-11-13 23:03:08 +01:00
Sven Czarnian
a8f2dba40e use the new interface to create the reports 2021-11-13 23:02:55 +01:00
Sven Czarnian
2b4840d2af fix a bug in the dynamics 2021-11-13 23:01:30 +01:00
Sven Czarnian
a3ef28ae39 rename the sequence receiver 2021-11-13 23:01:12 +01:00
Sven Czarnian
8b1292c95a rename the files and use the new interface to communicate with the backend 2021-11-13 22:59:04 +01:00
Sven Czarnian
50bb565229 build the new protocol files 2021-11-13 22:56:34 +01:00
Sven Czarnian
ab89237a4c extend the format 2021-10-11 21:56:13 +02:00
Sven Czarnian
02830454e7 fix the code to find the correct arrival route 2021-10-11 09:47:01 +02:00
Sven Czarnian
9f6f7c60d0 add the seconds to the report time 2021-10-11 09:46:47 +02:00
Sven Czarnian
9ce9503e7b update the submodule revision 2021-10-10 21:32:40 +02:00
Sven Czarnian
f4cef56164 calculate the required values 2021-10-10 21:32:29 +02:00
Sven Czarnian
818135e294 fix compiler warnings 2021-10-10 21:32:15 +02:00
Sven Czarnian
a934767ffb add the reported timestamp 2021-10-10 21:32:05 +02:00
Sven Czarnian
70a41bfc0e change the interface 2021-10-10 20:52:37 +02:00
Sven Czarnian
e0ab3ac68e fix compiler warnings 2021-10-10 20:52:20 +02:00
Sven Czarnian
7d69203486 updated protobuf to 3.18.1 2021-10-10 20:51:05 +02:00
Sven Czarnian
ea460b1aaf fix path errors 2021-10-10 20:50:20 +02:00
Sven Czarnian
8230538325 introduce the receiver of the scheduling information 2021-08-19 08:48:08 +02:00
Sven Czarnian
559ab86d60 remove unused headers 2021-08-19 08:47:50 +02:00
Sven Czarnian
2cf49c00fc introduce the aircraft scheduler to receive new messages from the backend
threading-concepts are used to avoid unnecessary main thread action
2021-08-19 08:47:39 +02:00
Sven Czarnian
c6677761dd check if the reporter is initialized 2021-08-19 08:21:53 +02:00
Sven Czarnian
eb05256df7 remove temporary filter 2021-08-19 08:21:29 +02:00
Sven Czarnian
221d00c23a add missing documentation 2021-08-19 08:21:20 +02:00
Sven Czarnian
b453d4d8e2 add a function to check if a block is initialized 2021-08-19 08:21:08 +02:00
Sven Czarnian
0fa3559fc4 add missing documentation 2021-08-19 08:20:54 +02:00
Sven Czarnian
9e15b72c54 convert the message in the reporter 2021-08-19 08:20:20 +02:00
Sven Czarnian
d06ec49b2e fix a deadlock or sporadic crash during shutdown 2021-08-19 08:01:22 +02:00
Sven Czarnian
25c6a1fa86 send aircraft reports 2021-08-17 17:31:28 +02:00
Sven Czarnian
7f4f967d83 build all new components 2021-08-17 17:31:08 +02:00
Sven Czarnian
8e36b094a5 introduce the aircraft reporter 2021-08-17 17:30:49 +02:00
Sven Czarnian
b1ab52d717 update the configuration 2021-08-17 17:30:10 +02:00
Sven Czarnian
627361035c define a context for all components 2021-08-17 17:25:59 +02:00
Sven Czarnian
c03b75fc89 update the interface 2021-08-17 17:25:46 +02:00
Sven Czarnian
804213317c add dependent libraries 2021-08-16 16:46:50 +02:00
Sven Czarnian
5a5825e1e4 fix imported target issues 2021-08-16 08:15:19 +02:00
Sven Czarnian
87ba0be57d extend the configuration 2021-08-16 08:08:41 +02:00
Sven Czarnian
db3d3bf48d update the communication version 2021-08-16 07:59:52 +02:00
Sven Czarnian
7b94e22121 add some more filters 2021-08-14 15:30:38 +02:00
Sven Czarnian
ed8fe7ef9a update the submodule 2021-08-14 15:16:01 +02:00
Sven Czarnian
d94e85b748 start with the interpretation of the position updates of aircrafts 2021-08-14 15:03:21 +02:00
Sven Czarnian
24f2ab45e3 generate and build the protobuf files and link against protobuf 2021-08-14 15:02:59 +02:00
Sven Czarnian
93c02cac91 add some spacings 2021-08-14 15:02:21 +02:00
Sven Czarnian
c350298fb5 update the submodule hash 2021-08-14 15:01:58 +02:00
Sven Czarnian
57a4aae223 add aman-com as a submodule 2021-08-14 14:52:48 +02:00
Sven Czarnian
8a3c2d3cb8 prepare the receive aircraft position updates 2021-08-14 10:37:08 +02:00
Sven Czarnian
c7ab3bd551 introduce the radar screen to receive aircraft information 2021-08-14 10:36:43 +02:00
Sven Czarnian
5da73b0b52 Merge branch 'feature/configuration' into 'develop'
Feature/configuration

See merge request nav/aman-es!2
2021-08-11 06:31:08 +00:00
Sven Czarnian
f188b85dcd load the configuration files 2021-08-11 08:29:38 +02:00
Sven Czarnian
fad3ca0de1 add a log message and find the path of the DLL to load the configurations later 2021-08-11 08:29:14 +02:00
Sven Czarnian
ac022f55fc add the configuration to the plugin 2021-08-11 08:27:51 +02:00
Sven Czarnian
272fb3aab0 add the communication file format 2021-08-10 08:42:53 +02:00
Sven Czarnian
3c4acffbb9 define the configuration to handle communication settings 2021-08-10 08:42:36 +02:00
Sven Czarnian
81edcdcea2 add a parser to read the UID per user 2021-08-10 08:42:00 +02:00
Sven Czarnian
163786698c define the base class to define different file formats 2021-08-10 08:40:53 +02:00
Sven Czarnian
5543d28bbf add GSL as a dependency 2021-08-10 08:39:45 +02:00
Sven Czarnian
4de1e2fc43 add an helper class to process strings 2021-08-10 08:39:34 +02:00
Sven Czarnian
557c781fa6 add GSL as a development dependency 2021-08-10 08:38:48 +02:00
Sven Czarnian
2a0dad3c7b update comments 2021-08-10 08:09:39 +02:00
Sven Czarnian
d594fd3252 remove symbols 2021-08-10 07:56:50 +02:00
Sven Czarnian
f0db197d7d remove GDI+ 2021-08-10 07:56:41 +02:00
Sven Czarnian
9696a66aab rename the namespace 2021-08-10 07:56:31 +02:00
Sven Czarnian
c7cf982671 remove unused variables 2021-08-10 07:20:14 +02:00
Sven Czarnian
a043314902 build the plugin 2021-08-09 22:28:49 +02:00
Sven Czarnian
40fa70186b add a missing include directory 2021-08-09 22:28:39 +02:00
Sven Czarnian
7addd2ab43 configure the build files 2021-08-09 22:28:26 +02:00
Sven Czarnian
328420c185 build the plugin 2021-08-09 22:28:09 +02:00
Sven Czarnian
18d285031c introduce the plugin definition 2021-08-09 22:28:00 +02:00
Sven Czarnian
0035c97c20 introduce resource files 2021-08-09 22:27:42 +02:00
Sven Czarnian
52870429a6 add missing templates for the automated file generation 2021-08-09 22:27:27 +02:00
Sven Czarnian
6caf47886a define the root-cmake-script 2021-08-09 22:07:28 +02:00
Sven Czarnian
27857b0369 introduce some required cmake scripts 2021-08-09 22:07:14 +02:00
Sven Czarnian
4cde0f9baf add external dependencies in a pre-built way to avoid incompatibilities 2021-08-09 22:07:00 +02:00
Sven Czarnian
5f702016ad add the license of the project 2021-08-09 22:06:03 +02:00
Sven Czarnian
539987868c update the documentation 2021-08-09 22:05:52 +02:00
801 changed files with 278507 additions and 1 deletions

4
.gitmodules vendored Normal file
View File

@@ -0,0 +1,4 @@
[submodule "src/com/protobuf"]
path = src/com/protobuf
url = git@git.vatsim-germany.org:nav/aman-com.git
branch = feature/protobuf

67
CMakeLists.txt Normal file
View File

@@ -0,0 +1,67 @@
# Author:
# Sven Czarnian <devel@svcz.de>
# Copyright:
# 2021 Sven Czarnian
# License:
# GPLv3
# Brief:
# Creates the AMAN-EuroScope solution
CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
# define the project information
PROJECT(ArrivalMANager LANGUAGES CXX VERSION "0.1.0")
# define the language parameters
SET_PROPERTY(GLOBAL PROPERTY USE_FOLDERS ON)
SET(CMAKE_CXX_STANDARD 20)
SET(CMAKE_CXX_STANDARD_REQUIRED ON)
SET(CMAKE_CXX_EXTENSIONS OFF)
SET(CMAKE_INTERPROCEDURAL_OPTIMIZATION_RELEASE TRUE)
# adapt compiler flags based on used compiler
IF (MSVC)
IF (CMAKE_CXX_FLAGS MATCHES "/W[0-4]")
STRING(REGEX REPLACE "/W[0-4]" "/W4" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
ELSE ()
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4")
ENDIF ()
IF (NOT CMAKE_CXX_FLAGS MATCHES "/MP")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP")
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /MP")
ENDIF ()
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /sdl /permissive- /DNOMINMAX")
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /sdl /permissive- /DNOMINMAX")
SET(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} /MANIFESTUAC:NO")
ADD_DEFINITIONS(/D_USRDLL /D_CRT_SECURE_NO_WARNINGS)
ENDIF ()
CONFIGURE_FILE(
${CMAKE_SOURCE_DIR}/version.h.in
${CMAKE_BINARY_DIR}/include/version.h
)
CONFIGURE_FILE(
${CMAKE_SOURCE_DIR}/res/ArrivalMANager.rc.in
${CMAKE_BINARY_DIR}/ArrivalMANager.rc
)
# define include directories
INCLUDE_DIRECTORIES(
${CMAKE_SOURCE_DIR}
${CMAKE_SOURCE_DIR}/include
${CMAKE_BINARY_DIR}/include
)
INCLUDE(cmake/3rdParty.cmake)
INCLUDE(cmake/Protobuf.cmake)
INCLUDE(cmake/FindEuroScope.cmake)
# register all cmake helper to find required modules and find 3rd-party components
SET(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${CMAKE_SOURCE_DIR}/cmake")
FIND_PACKAGE(EuroScope REQUIRED)
IF(NOT EuroScope_FOUND)
MESSAGE(FATAL_ERROR "Unablet to build without EuroScope and the EuroScope-SDK")
ENDIF()
ADD_SUBDIRECTORY(src)

674
LICENSE Normal file
View File

@@ -0,0 +1,674 @@
GNU GENERAL PUBLIC LICENSE
Version 3, 29 June 2007
Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
Preamble
The GNU General Public License is a free, copyleft license for
software and other kinds of works.
The licenses for most software and other practical works are designed
to take away your freedom to share and change the works. By contrast,
the GNU General Public License is intended to guarantee your freedom to
share and change all versions of a program--to make sure it remains free
software for all its users. We, the Free Software Foundation, use the
GNU General Public License for most of our software; it applies also to
any other work released this way by its authors. You can apply it to
your programs, too.
When we speak of free software, we are referring to freedom, not
price. Our General Public Licenses are designed to make sure that you
have the freedom to distribute copies of free software (and charge for
them if you wish), that you receive source code or can get it if you
want it, that you can change the software or use pieces of it in new
free programs, and that you know you can do these things.
To protect your rights, we need to prevent others from denying you
these rights or asking you to surrender the rights. Therefore, you have
certain responsibilities if you distribute copies of the software, or if
you modify it: responsibilities to respect the freedom of others.
For example, if you distribute copies of such a program, whether
gratis or for a fee, you must pass on to the recipients the same
freedoms that you received. You must make sure that they, too, receive
or can get the source code. And you must show them these terms so they
know their rights.
Developers that use the GNU GPL protect your rights with two steps:
(1) assert copyright on the software, and (2) offer you this License
giving you legal permission to copy, distribute and/or modify it.
For the developers' and authors' protection, the GPL clearly explains
that there is no warranty for this free software. For both users' and
authors' sake, the GPL requires that modified versions be marked as
changed, so that their problems will not be attributed erroneously to
authors of previous versions.
Some devices are designed to deny users access to install or run
modified versions of the software inside them, although the manufacturer
can do so. This is fundamentally incompatible with the aim of
protecting users' freedom to change the software. The systematic
pattern of such abuse occurs in the area of products for individuals to
use, which is precisely where it is most unacceptable. Therefore, we
have designed this version of the GPL to prohibit the practice for those
products. If such problems arise substantially in other domains, we
stand ready to extend this provision to those domains in future versions
of the GPL, as needed to protect the freedom of users.
Finally, every program is threatened constantly by software patents.
States should not allow patents to restrict development and use of
software on general-purpose computers, but in those that do, we wish to
avoid the special danger that patents applied to a free program could
make it effectively proprietary. To prevent this, the GPL assures that
patents cannot be used to render the program non-free.
The precise terms and conditions for copying, distribution and
modification follow.
TERMS AND CONDITIONS
0. Definitions.
"This License" refers to version 3 of the GNU General Public License.
"Copyright" also means copyright-like laws that apply to other kinds of
works, such as semiconductor masks.
"The Program" refers to any copyrightable work licensed under this
License. Each licensee is addressed as "you". "Licensees" and
"recipients" may be individuals or organizations.
To "modify" a work means to copy from or adapt all or part of the work
in a fashion requiring copyright permission, other than the making of an
exact copy. The resulting work is called a "modified version" of the
earlier work or a work "based on" the earlier work.
A "covered work" means either the unmodified Program or a work based
on the Program.
To "propagate" a work means to do anything with it that, without
permission, would make you directly or secondarily liable for
infringement under applicable copyright law, except executing it on a
computer or modifying a private copy. Propagation includes copying,
distribution (with or without modification), making available to the
public, and in some countries other activities as well.
To "convey" a work means any kind of propagation that enables other
parties to make or receive copies. Mere interaction with a user through
a computer network, with no transfer of a copy, is not conveying.
An interactive user interface displays "Appropriate Legal Notices"
to the extent that it includes a convenient and prominently visible
feature that (1) displays an appropriate copyright notice, and (2)
tells the user that there is no warranty for the work (except to the
extent that warranties are provided), that licensees may convey the
work under this License, and how to view a copy of this License. If
the interface presents a list of user commands or options, such as a
menu, a prominent item in the list meets this criterion.
1. Source Code.
The "source code" for a work means the preferred form of the work
for making modifications to it. "Object code" means any non-source
form of a work.
A "Standard Interface" means an interface that either is an official
standard defined by a recognized standards body, or, in the case of
interfaces specified for a particular programming language, one that
is widely used among developers working in that language.
The "System Libraries" of an executable work include anything, other
than the work as a whole, that (a) is included in the normal form of
packaging a Major Component, but which is not part of that Major
Component, and (b) serves only to enable use of the work with that
Major Component, or to implement a Standard Interface for which an
implementation is available to the public in source code form. A
"Major Component", in this context, means a major essential component
(kernel, window system, and so on) of the specific operating system
(if any) on which the executable work runs, or a compiler used to
produce the work, or an object code interpreter used to run it.
The "Corresponding Source" for a work in object code form means all
the source code needed to generate, install, and (for an executable
work) run the object code and to modify the work, including scripts to
control those activities. However, it does not include the work's
System Libraries, or general-purpose tools or generally available free
programs which are used unmodified in performing those activities but
which are not part of the work. For example, Corresponding Source
includes interface definition files associated with source files for
the work, and the source code for shared libraries and dynamically
linked subprograms that the work is specifically designed to require,
such as by intimate data communication or control flow between those
subprograms and other parts of the work.
The Corresponding Source need not include anything that users
can regenerate automatically from other parts of the Corresponding
Source.
The Corresponding Source for a work in source code form is that
same work.
2. Basic Permissions.
All rights granted under this License are granted for the term of
copyright on the Program, and are irrevocable provided the stated
conditions are met. This License explicitly affirms your unlimited
permission to run the unmodified Program. The output from running a
covered work is covered by this License only if the output, given its
content, constitutes a covered work. This License acknowledges your
rights of fair use or other equivalent, as provided by copyright law.
You may make, run and propagate covered works that you do not
convey, without conditions so long as your license otherwise remains
in force. You may convey covered works to others for the sole purpose
of having them make modifications exclusively for you, or provide you
with facilities for running those works, provided that you comply with
the terms of this License in conveying all material for which you do
not control copyright. Those thus making or running the covered works
for you must do so exclusively on your behalf, under your direction
and control, on terms that prohibit them from making any copies of
your copyrighted material outside their relationship with you.
Conveying under any other circumstances is permitted solely under
the conditions stated below. Sublicensing is not allowed; section 10
makes it unnecessary.
3. Protecting Users' Legal Rights From Anti-Circumvention Law.
No covered work shall be deemed part of an effective technological
measure under any applicable law fulfilling obligations under article
11 of the WIPO copyright treaty adopted on 20 December 1996, or
similar laws prohibiting or restricting circumvention of such
measures.
When you convey a covered work, you waive any legal power to forbid
circumvention of technological measures to the extent such circumvention
is effected by exercising rights under this License with respect to
the covered work, and you disclaim any intention to limit operation or
modification of the work as a means of enforcing, against the work's
users, your or third parties' legal rights to forbid circumvention of
technological measures.
4. Conveying Verbatim Copies.
You may convey verbatim copies of the Program's source code as you
receive it, in any medium, provided that you conspicuously and
appropriately publish on each copy an appropriate copyright notice;
keep intact all notices stating that this License and any
non-permissive terms added in accord with section 7 apply to the code;
keep intact all notices of the absence of any warranty; and give all
recipients a copy of this License along with the Program.
You may charge any price or no price for each copy that you convey,
and you may offer support or warranty protection for a fee.
5. Conveying Modified Source Versions.
You may convey a work based on the Program, or the modifications to
produce it from the Program, in the form of source code under the
terms of section 4, provided that you also meet all of these conditions:
a) The work must carry prominent notices stating that you modified
it, and giving a relevant date.
b) The work must carry prominent notices stating that it is
released under this License and any conditions added under section
7. This requirement modifies the requirement in section 4 to
"keep intact all notices".
c) You must license the entire work, as a whole, under this
License to anyone who comes into possession of a copy. This
License will therefore apply, along with any applicable section 7
additional terms, to the whole of the work, and all its parts,
regardless of how they are packaged. This License gives no
permission to license the work in any other way, but it does not
invalidate such permission if you have separately received it.
d) If the work has interactive user interfaces, each must display
Appropriate Legal Notices; however, if the Program has interactive
interfaces that do not display Appropriate Legal Notices, your
work need not make them do so.
A compilation of a covered work with other separate and independent
works, which are not by their nature extensions of the covered work,
and which are not combined with it such as to form a larger program,
in or on a volume of a storage or distribution medium, is called an
"aggregate" if the compilation and its resulting copyright are not
used to limit the access or legal rights of the compilation's users
beyond what the individual works permit. Inclusion of a covered work
in an aggregate does not cause this License to apply to the other
parts of the aggregate.
6. Conveying Non-Source Forms.
You may convey a covered work in object code form under the terms
of sections 4 and 5, provided that you also convey the
machine-readable Corresponding Source under the terms of this License,
in one of these ways:
a) Convey the object code in, or embodied in, a physical product
(including a physical distribution medium), accompanied by the
Corresponding Source fixed on a durable physical medium
customarily used for software interchange.
b) Convey the object code in, or embodied in, a physical product
(including a physical distribution medium), accompanied by a
written offer, valid for at least three years and valid for as
long as you offer spare parts or customer support for that product
model, to give anyone who possesses the object code either (1) a
copy of the Corresponding Source for all the software in the
product that is covered by this License, on a durable physical
medium customarily used for software interchange, for a price no
more than your reasonable cost of physically performing this
conveying of source, or (2) access to copy the
Corresponding Source from a network server at no charge.
c) Convey individual copies of the object code with a copy of the
written offer to provide the Corresponding Source. This
alternative is allowed only occasionally and noncommercially, and
only if you received the object code with such an offer, in accord
with subsection 6b.
d) Convey the object code by offering access from a designated
place (gratis or for a charge), and offer equivalent access to the
Corresponding Source in the same way through the same place at no
further charge. You need not require recipients to copy the
Corresponding Source along with the object code. If the place to
copy the object code is a network server, the Corresponding Source
may be on a different server (operated by you or a third party)
that supports equivalent copying facilities, provided you maintain
clear directions next to the object code saying where to find the
Corresponding Source. Regardless of what server hosts the
Corresponding Source, you remain obligated to ensure that it is
available for as long as needed to satisfy these requirements.
e) Convey the object code using peer-to-peer transmission, provided
you inform other peers where the object code and Corresponding
Source of the work are being offered to the general public at no
charge under subsection 6d.
A separable portion of the object code, whose source code is excluded
from the Corresponding Source as a System Library, need not be
included in conveying the object code work.
A "User Product" is either (1) a "consumer product", which means any
tangible personal property which is normally used for personal, family,
or household purposes, or (2) anything designed or sold for incorporation
into a dwelling. In determining whether a product is a consumer product,
doubtful cases shall be resolved in favor of coverage. For a particular
product received by a particular user, "normally used" refers to a
typical or common use of that class of product, regardless of the status
of the particular user or of the way in which the particular user
actually uses, or expects or is expected to use, the product. A product
is a consumer product regardless of whether the product has substantial
commercial, industrial or non-consumer uses, unless such uses represent
the only significant mode of use of the product.
"Installation Information" for a User Product means any methods,
procedures, authorization keys, or other information required to install
and execute modified versions of a covered work in that User Product from
a modified version of its Corresponding Source. The information must
suffice to ensure that the continued functioning of the modified object
code is in no case prevented or interfered with solely because
modification has been made.
If you convey an object code work under this section in, or with, or
specifically for use in, a User Product, and the conveying occurs as
part of a transaction in which the right of possession and use of the
User Product is transferred to the recipient in perpetuity or for a
fixed term (regardless of how the transaction is characterized), the
Corresponding Source conveyed under this section must be accompanied
by the Installation Information. But this requirement does not apply
if neither you nor any third party retains the ability to install
modified object code on the User Product (for example, the work has
been installed in ROM).
The requirement to provide Installation Information does not include a
requirement to continue to provide support service, warranty, or updates
for a work that has been modified or installed by the recipient, or for
the User Product in which it has been modified or installed. Access to a
network may be denied when the modification itself materially and
adversely affects the operation of the network or violates the rules and
protocols for communication across the network.
Corresponding Source conveyed, and Installation Information provided,
in accord with this section must be in a format that is publicly
documented (and with an implementation available to the public in
source code form), and must require no special password or key for
unpacking, reading or copying.
7. Additional Terms.
"Additional permissions" are terms that supplement the terms of this
License by making exceptions from one or more of its conditions.
Additional permissions that are applicable to the entire Program shall
be treated as though they were included in this License, to the extent
that they are valid under applicable law. If additional permissions
apply only to part of the Program, that part may be used separately
under those permissions, but the entire Program remains governed by
this License without regard to the additional permissions.
When you convey a copy of a covered work, you may at your option
remove any additional permissions from that copy, or from any part of
it. (Additional permissions may be written to require their own
removal in certain cases when you modify the work.) You may place
additional permissions on material, added by you to a covered work,
for which you have or can give appropriate copyright permission.
Notwithstanding any other provision of this License, for material you
add to a covered work, you may (if authorized by the copyright holders of
that material) supplement the terms of this License with terms:
a) Disclaiming warranty or limiting liability differently from the
terms of sections 15 and 16 of this License; or
b) Requiring preservation of specified reasonable legal notices or
author attributions in that material or in the Appropriate Legal
Notices displayed by works containing it; or
c) Prohibiting misrepresentation of the origin of that material, or
requiring that modified versions of such material be marked in
reasonable ways as different from the original version; or
d) Limiting the use for publicity purposes of names of licensors or
authors of the material; or
e) Declining to grant rights under trademark law for use of some
trade names, trademarks, or service marks; or
f) Requiring indemnification of licensors and authors of that
material by anyone who conveys the material (or modified versions of
it) with contractual assumptions of liability to the recipient, for
any liability that these contractual assumptions directly impose on
those licensors and authors.
All other non-permissive additional terms are considered "further
restrictions" within the meaning of section 10. If the Program as you
received it, or any part of it, contains a notice stating that it is
governed by this License along with a term that is a further
restriction, you may remove that term. If a license document contains
a further restriction but permits relicensing or conveying under this
License, you may add to a covered work material governed by the terms
of that license document, provided that the further restriction does
not survive such relicensing or conveying.
If you add terms to a covered work in accord with this section, you
must place, in the relevant source files, a statement of the
additional terms that apply to those files, or a notice indicating
where to find the applicable terms.
Additional terms, permissive or non-permissive, may be stated in the
form of a separately written license, or stated as exceptions;
the above requirements apply either way.
8. Termination.
You may not propagate or modify a covered work except as expressly
provided under this License. Any attempt otherwise to propagate or
modify it is void, and will automatically terminate your rights under
this License (including any patent licenses granted under the third
paragraph of section 11).
However, if you cease all violation of this License, then your
license from a particular copyright holder is reinstated (a)
provisionally, unless and until the copyright holder explicitly and
finally terminates your license, and (b) permanently, if the copyright
holder fails to notify you of the violation by some reasonable means
prior to 60 days after the cessation.
Moreover, your license from a particular copyright holder is
reinstated permanently if the copyright holder notifies you of the
violation by some reasonable means, this is the first time you have
received notice of violation of this License (for any work) from that
copyright holder, and you cure the violation prior to 30 days after
your receipt of the notice.
Termination of your rights under this section does not terminate the
licenses of parties who have received copies or rights from you under
this License. If your rights have been terminated and not permanently
reinstated, you do not qualify to receive new licenses for the same
material under section 10.
9. Acceptance Not Required for Having Copies.
You are not required to accept this License in order to receive or
run a copy of the Program. Ancillary propagation of a covered work
occurring solely as a consequence of using peer-to-peer transmission
to receive a copy likewise does not require acceptance. However,
nothing other than this License grants you permission to propagate or
modify any covered work. These actions infringe copyright if you do
not accept this License. Therefore, by modifying or propagating a
covered work, you indicate your acceptance of this License to do so.
10. Automatic Licensing of Downstream Recipients.
Each time you convey a covered work, the recipient automatically
receives a license from the original licensors, to run, modify and
propagate that work, subject to this License. You are not responsible
for enforcing compliance by third parties with this License.
An "entity transaction" is a transaction transferring control of an
organization, or substantially all assets of one, or subdividing an
organization, or merging organizations. If propagation of a covered
work results from an entity transaction, each party to that
transaction who receives a copy of the work also receives whatever
licenses to the work the party's predecessor in interest had or could
give under the previous paragraph, plus a right to possession of the
Corresponding Source of the work from the predecessor in interest, if
the predecessor has it or can get it with reasonable efforts.
You may not impose any further restrictions on the exercise of the
rights granted or affirmed under this License. For example, you may
not impose a license fee, royalty, or other charge for exercise of
rights granted under this License, and you may not initiate litigation
(including a cross-claim or counterclaim in a lawsuit) alleging that
any patent claim is infringed by making, using, selling, offering for
sale, or importing the Program or any portion of it.
11. Patents.
A "contributor" is a copyright holder who authorizes use under this
License of the Program or a work on which the Program is based. The
work thus licensed is called the contributor's "contributor version".
A contributor's "essential patent claims" are all patent claims
owned or controlled by the contributor, whether already acquired or
hereafter acquired, that would be infringed by some manner, permitted
by this License, of making, using, or selling its contributor version,
but do not include claims that would be infringed only as a
consequence of further modification of the contributor version. For
purposes of this definition, "control" includes the right to grant
patent sublicenses in a manner consistent with the requirements of
this License.
Each contributor grants you a non-exclusive, worldwide, royalty-free
patent license under the contributor's essential patent claims, to
make, use, sell, offer for sale, import and otherwise run, modify and
propagate the contents of its contributor version.
In the following three paragraphs, a "patent license" is any express
agreement or commitment, however denominated, not to enforce a patent
(such as an express permission to practice a patent or covenant not to
sue for patent infringement). To "grant" such a patent license to a
party means to make such an agreement or commitment not to enforce a
patent against the party.
If you convey a covered work, knowingly relying on a patent license,
and the Corresponding Source of the work is not available for anyone
to copy, free of charge and under the terms of this License, through a
publicly available network server or other readily accessible means,
then you must either (1) cause the Corresponding Source to be so
available, or (2) arrange to deprive yourself of the benefit of the
patent license for this particular work, or (3) arrange, in a manner
consistent with the requirements of this License, to extend the patent
license to downstream recipients. "Knowingly relying" means you have
actual knowledge that, but for the patent license, your conveying the
covered work in a country, or your recipient's use of the covered work
in a country, would infringe one or more identifiable patents in that
country that you have reason to believe are valid.
If, pursuant to or in connection with a single transaction or
arrangement, you convey, or propagate by procuring conveyance of, a
covered work, and grant a patent license to some of the parties
receiving the covered work authorizing them to use, propagate, modify
or convey a specific copy of the covered work, then the patent license
you grant is automatically extended to all recipients of the covered
work and works based on it.
A patent license is "discriminatory" if it does not include within
the scope of its coverage, prohibits the exercise of, or is
conditioned on the non-exercise of one or more of the rights that are
specifically granted under this License. You may not convey a covered
work if you are a party to an arrangement with a third party that is
in the business of distributing software, under which you make payment
to the third party based on the extent of your activity of conveying
the work, and under which the third party grants, to any of the
parties who would receive the covered work from you, a discriminatory
patent license (a) in connection with copies of the covered work
conveyed by you (or copies made from those copies), or (b) primarily
for and in connection with specific products or compilations that
contain the covered work, unless you entered into that arrangement,
or that patent license was granted, prior to 28 March 2007.
Nothing in this License shall be construed as excluding or limiting
any implied license or other defenses to infringement that may
otherwise be available to you under applicable patent law.
12. No Surrender of Others' Freedom.
If conditions are imposed on you (whether by court order, agreement or
otherwise) that contradict the conditions of this License, they do not
excuse you from the conditions of this License. If you cannot convey a
covered work so as to satisfy simultaneously your obligations under this
License and any other pertinent obligations, then as a consequence you may
not convey it at all. For example, if you agree to terms that obligate you
to collect a royalty for further conveying from those to whom you convey
the Program, the only way you could satisfy both those terms and this
License would be to refrain entirely from conveying the Program.
13. Use with the GNU Affero General Public License.
Notwithstanding any other provision of this License, you have
permission to link or combine any covered work with a work licensed
under version 3 of the GNU Affero General Public License into a single
combined work, and to convey the resulting work. The terms of this
License will continue to apply to the part which is the covered work,
but the special requirements of the GNU Affero General Public License,
section 13, concerning interaction through a network will apply to the
combination as such.
14. Revised Versions of this License.
The Free Software Foundation may publish revised and/or new versions of
the GNU General Public License from time to time. Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.
Each version is given a distinguishing version number. If the
Program specifies that a certain numbered version of the GNU General
Public License "or any later version" applies to it, you have the
option of following the terms and conditions either of that numbered
version or of any later version published by the Free Software
Foundation. If the Program does not specify a version number of the
GNU General Public License, you may choose any version ever published
by the Free Software Foundation.
If the Program specifies that a proxy can decide which future
versions of the GNU General Public License can be used, that proxy's
public statement of acceptance of a version permanently authorizes you
to choose that version for the Program.
Later license versions may give you additional or different
permissions. However, no additional obligations are imposed on any
author or copyright holder as a result of your choosing to follow a
later version.
15. Disclaimer of Warranty.
THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
16. Limitation of Liability.
IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
SUCH DAMAGES.
17. Interpretation of Sections 15 and 16.
If the disclaimer of warranty and limitation of liability provided
above cannot be given local legal effect according to their terms,
reviewing courts shall apply local law that most closely approximates
an absolute waiver of all civil liability in connection with the
Program, unless a warranty or assumption of liability accompanies a
copy of the Program in return for a fee.
END OF TERMS AND CONDITIONS
How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
free software which everyone can redistribute and change under these terms.
To do so, attach the following notices to the program. It is safest
to attach them to the start of each source file to most effectively
state the exclusion of warranty; and each file should have at least
the "copyright" line and a pointer to where the full notice is found.
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
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 3 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, see <https://www.gnu.org/licenses/>.
Also add information on how to contact you by electronic and paper mail.
If the program does terminal interaction, make it output a short
notice like this when it starts in an interactive mode:
<program> Copyright (C) <year> <name of author>
This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
The hypothetical commands `show w' and `show c' should show the appropriate
parts of the General Public License. Of course, your program's commands
might be different; for a GUI interface, you would use an "about box".
You should also get your employer (if you work as a programmer) or school,
if any, to sign a "copyright disclaimer" for the program, if necessary.
For more information on this, and how to apply and follow the GNU GPL, see
<https://www.gnu.org/licenses/>.
The GNU General Public License does not permit incorporating your program
into proprietary programs. If your program is a subroutine library, you
may consider it more useful to permit linking proprietary applications with
the library. If this is what you want to do, use the GNU Lesser General
Public License instead of this License. But first, please read
<https://www.gnu.org/licenses/why-not-lgpl.html>.

View File

@@ -1,2 +1,30 @@
# aman-es
# Arrival MANanager (AMAN)
## System description
AMAN is splitted up into four different components.
* aman-com defines the diffent message types
* aman-es implements an EuroScope plugin to communicate with [aman-sys](https://git.vatsim-germany.org/nav/aman-sys)
* aman-sys implements the backend system to plan an optimal arrival sequence for the different airports
* aman-web implements a web-interface to configure [aman-sys](https://git.vatsim-germany.org/nav/aman-sys) and visualize sequences
## Component description
AMAN uses [Protocol Buffers](https://developers.google.com/protocol-buffers)
for the message serialization and message definition between the EuroScope instance and the AMAN backend.
Additionally is [ZeroMQ](https://zeromq.org/) used for the communication abstraction layer.
This component extracts all relevant information out of the VATSIM network and sends the information to the server.
Additionally is a unique identifier used to verify connections to the backend.
According to data protection is the ZeroMQ-based network encryption used.
Every controller needs his own unique identifier.
# Additional libraries
* [ZeroMQ](https://github.com/zeromq) - GNU GPLv3
* [Protocol Buffers](https://github.com/protocolbuffers/protobuf) - BSD-3
# License
AMAN is released under the [GNU General Public License v3](LICENSE)

8
cmake/3rdParty.cmake Normal file
View File

@@ -0,0 +1,8 @@
# Author:
# Sven Czarnian <devel@svcz.de>
# License:
# GPLv3
# Brief:
# Creates the 3rd-party targets
INCLUDE(${CMAKE_SOURCE_DIR}/cmake/3rdPartyTargets.cmake)

100
cmake/3rdPartyTargets.cmake Normal file
View File

@@ -0,0 +1,100 @@
# Author:
# Sven Czarnian <devel@svcz.de>
# License:
# GPLv3
# Brief:
# Creates the 3rd-party import targets
# include the external project library
INCLUDE(ExternalProject)
# define the import target of libcurl
ADD_LIBRARY(libcurl STATIC IMPORTED)
IF (MSVC)
SET_TARGET_PROPERTIES(libcurl PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/libcurl-d.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/libcurl.lib"
)
TARGET_LINK_LIBRARIES(libcurl INTERFACE Ws2_32 Iphlpapi Crypt32)
TARGET_INCLUDE_DIRECTORIES(libcurl INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()
# define the import target of GeographicLib
ADD_LIBRARY(GeographicLib STATIC IMPORTED)
IF (MSVC)
SET_TARGET_PROPERTIES(GeographicLib PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/Geographic_d.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/Geographic.lib"
)
TARGET_INCLUDE_DIRECTORIES(GeographicLib INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()
# define the import target of libsodium
ADD_LIBRARY(libsodium STATIC IMPORTED)
IF (MSVC)
SET_TARGET_PROPERTIES(libsodium PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/libsodiumd.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/libsodium.lib"
)
TARGET_LINK_LIBRARIES(libsodium INTERFACE Ws2_32 Iphlpapi)
TARGET_INCLUDE_DIRECTORIES(libsodium INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()
# define the import target of libcurl
ADD_LIBRARY(jsoncpp STATIC IMPORTED)
IF (MSVC)
SET_TARGET_PROPERTIES(jsoncpp PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/jsoncpp_d.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/jsoncpp.lib"
)
TARGET_INCLUDE_DIRECTORIES(jsoncpp INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()
# define the import target of libzmq
ADD_LIBRARY(libzmq STATIC IMPORTED)
ADD_DEPENDENCIES(libzmq libsodium)
IF (MSVC)
SET_TARGET_PROPERTIES(libzmq PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/libzmqd.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/libzmq.lib"
)
TARGET_INCLUDE_DIRECTORIES(libzmq INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
TARGET_LINK_LIBRARIES(libzmq INTERFACE libsodium)
TARGET_COMPILE_OPTIONS(libzmq INTERFACE /DZMQ_STATIC)
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()
# define the import target of cppzmq
ADD_LIBRARY(cppzmq INTERFACE)
TARGET_INCLUDE_DIRECTORIES(cppzmq INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
TARGET_LINK_LIBRARIES(cppzmq INTERFACE libzmq)
ADD_DEPENDENCIES(cppzmq libzmq)
# define the import target of GSL
ADD_LIBRARY(GSL INTERFACE)
TARGET_INCLUDE_DIRECTORIES(GSL INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
# define the import target of Eigen
ADD_LIBRARY(Eigen INTERFACE)
TARGET_INCLUDE_DIRECTORIES(Eigen INTERFACE "${CMAKE_SOURCE_DIR}/external/include/eigen3")
# define the import target of protobuf
ADD_LIBRARY(protobuf STATIC IMPORTED)
IF (MSVC)
SET_TARGET_PROPERTIES(protobuf PROPERTIES
IMPORTED_LOCATION_DEBUG "${CMAKE_SOURCE_DIR}/external/lib/libprotobufd.lib"
IMPORTED_LOCATION_RELEASE "${CMAKE_SOURCE_DIR}/external/lib/libprotobuf.lib"
)
TARGET_INCLUDE_DIRECTORIES(protobuf INTERFACE "${CMAKE_SOURCE_DIR}/external/include")
ELSE ()
MESSAGE(FATAL_ERROR "Unsupported compiler")
ENDIF ()

49
cmake/FindEuroScope.cmake Normal file
View File

@@ -0,0 +1,49 @@
# Author:
# Sven Czarnian <devel@svcz.de>
# License:
# LGPLv3
# Brief:
# Finds the EuroScope headers and libraries
# A target EuroScope will be created and the EuroScope_FOUND flag will be set
IF(NOT TARGET EuroScope)
IF(NOT EuroScope_DIR)
MESSAGE(FATAL_ERROR "Please set EuroScope_DIR")
SET(EuroScope_DIR "EuroScope_DIR-NOTFOUND" CACHE PATH PARENT_SCOPE)
ENDIF()
FIND_FILE(EuroScope_EXECUTABLE
NAMES
EuroScope.exe
PATHS
${EuroScope_DIR}
)
FIND_FILE(EuroScope_LIBRARY
NAMES
EuroScopePlugInDll.lib
PATHS
${EuroScope_DIR}/PlugInEnvironment
)
FIND_PATH(EuroScope_INCLUDE_DIR
NAMES
EuroScopePlugIn.h
PATHS
${EuroScope_DIR}/PlugInEnvironment
)
IF(NOT ${EuroScope_EXECUTABLE} STREQUAL "EuroScope_EXECUTABLE-NOTFOUND" AND
NOT ${EuroScope_LIBRARY} STREQUAL "EuroScope_LIBRARY-NOTFOUND" AND
NOT ${EuroScope_INCLUDE_DIR} STREQUAL "EuroScope_INCLUDE_DIR-NOTFOUND")
MESSAGE(STATUS "Found EuroScope-library:")
MESSAGE(STATUS " ${EuroScope_LIBRARY}")
MESSAGE(STATUS "Found EuroScope-headers:")
MESSAGE(STATUS " ${EuroScope_INCLUDE_DIR}")
ADD_LIBRARY(EuroScope INTERFACE IMPORTED GLOBAL)
TARGET_LINK_LIBRARIES(EuroScope INTERFACE ${EuroScope_LIBRARY})
TARGET_INCLUDE_DIRECTORIES(EuroScope INTERFACE ${EuroScope_INCLUDE_DIR})
SET(EuroScope_FOUND ON)
ENDIF()
ELSE()
MESSAGE(STATUS "EuroScope is already included.")
ENDIF()

51
cmake/Protobuf.cmake Normal file
View File

@@ -0,0 +1,51 @@
# Author:
# Sven Czarnian <devel@svcz.de>
# License:
# Closed Source
# Brief:
# Defines the protobuf functions
# Brief:
# Proto-files are compiled into C++ files
# Parameters:
# PROTO_FILES - The proto-files with the message description
# SOURCE_FILES - Contains the filenames and paths of the generated files
FUNCTION(ProtobufCompile PROTO_FILES SOURCE_FILES)
SET(GENERATED_FILES "")
FOREACH (PROTO ${PROTO_FILES})
# get the relevant information to configure the protoc-run
GET_FILENAME_COMPONENT(FILENAME ${PROTO} NAME_WLE)
GET_FILENAME_COMPONENT(DIRECTORY ${PROTO} DIRECTORY)
# define the output files
SET(CPP_FILE ${CMAKE_CURRENT_BINARY_DIR}/protobuf/${FILENAME}.pb.cc)
SET(HPP_FILE ${CMAKE_CURRENT_BINARY_DIR}/protobuf/${FILENAME}.pb.h)
# create the protoc-directory
FILE(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/protobuf)
# define the protoc-command
ADD_CUSTOM_COMMAND(
OUTPUT ${CPP_FILE} ${HPP_FILE}
DEPENDS protobuf
COMMAND ${CMAKE_SOURCE_DIR}/external/bin/protoc.exe
ARGS -I=${DIRECTORY} --cpp_out=${CMAKE_CURRENT_BINARY_DIR}/protobuf ${PROTO}
WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}/external/bin"
COMMENT "Creating C++-sources for ${PROTO}"
)
# disable warnings
IF (MSVC)
SET_SOURCE_FILES_PROPERTIES(${CPP_FILE} PROPERTIES COMPILE_FLAGS "/wd4127 /wd5054 /wd4125 /wd4267")
SET_SOURCE_FILES_PROPERTIES(${HPP_FILE} PROPERTIES COMPILE_FLAGS "/wd4127 /wd5054 /wd4125 /wd4267")
ENDIF ()
# add the generated files
LIST(APPEND GENERATED_FILES ${CPP_FILE})
LIST(APPEND GENERATED_FILES ${HPP_FILE})
ENDFOREACH ()
# set the output variables
SET(${SOURCE_FILES} ${GENERATED_FILES} PARENT_SCOPE)
ENDFUNCTION()

BIN
external/bin/protoc.exe vendored Normal file

Binary file not shown.

View File

@@ -0,0 +1,198 @@
/**
* \file Accumulator.hpp
* \brief Header for GeographicLib::Accumulator class
*
* Copyright (c) Charles Karney (2010-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP)
#define GEOGRAPHICLIB_ACCUMULATOR_HPP 1
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief An accumulator for sums
*
* This allows many numbers of floating point type \e T to be added together
* with twice the normal precision. Thus if \e T is double, the effective
* precision of the sum is 106 bits or about 32 decimal places.
*
* The implementation follows J. R. Shewchuk,
* <a href="https://doi.org/10.1007/PL00009321"> Adaptive Precision
* Floating-Point Arithmetic and Fast Robust Geometric Predicates</a>,
* Discrete & Computational Geometry 18(3) 305--363 (1997).
*
* Approximate timings (summing a vector<double>)
* - double: 2ns
* - Accumulator<double>: 23ns
*
* In the documentation of the member functions, \e sum stands for the value
* currently held in the accumulator.
*
* Example of use:
* \include example-Accumulator.cpp
**********************************************************************/
template<typename T = Math::real>
class GEOGRAPHICLIB_EXPORT Accumulator {
private:
// _s + _t accumulators for the sum.
T _s, _t;
// Same as Math::sum, but requires abs(u) >= abs(v). This isn't currently
// used.
static T fastsum(T u, T v, T& t) {
GEOGRAPHICLIB_VOLATILE T s = u + v;
GEOGRAPHICLIB_VOLATILE T vp = s - u;
t = v - vp;
return s;
}
void Add(T y) {
// Here's Shewchuk's solution...
T u; // hold exact sum as [s, t, u]
// Accumulate starting at least significant end
y = Math::sum(y, _t, u);
_s = Math::sum(y, _s, _t);
// Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u)
// exactly with s, t, u non-adjacent and in decreasing order (except for
// possible zeros). The following code tries to normalize the result.
// Ideally, we want _s = round(s+t+u) and _u = round(s+t+u - _s). The
// following does an approximate job (and maintains the decreasing
// non-adjacent property). Here are two "failures" using 3-bit floats:
//
// Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
// [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
//
// Case 2: _s+_t is not as close to s+t+u as it shold be
// [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
// should be [80, -7] = 73 (exact)
//
// "Fixing" these problems is probably not worth the expense. The
// representation inevitably leads to small errors in the accumulated
// values. The additional errors illustrated here amount to 1 ulp of the
// less significant word during each addition to the Accumulator and an
// additional possible error of 1 ulp in the reported sum.
//
// Incidentally, the "ideal" representation described above is not
// canonical, because _s = round(_s + _t) may not be true. For example,
// with 3-bit floats:
//
// [128, 16] + 1 -> [160, -16] -- 160 = round(145).
// But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
//
if (_s == 0) // This implies t == 0,
_s = u; // so result is u
else
_t += u; // otherwise just accumulate u to t.
}
T Sum(T y) const {
Accumulator a(*this);
a.Add(y);
return a._s;
}
public:
/**
* Construct from a \e T. This is not declared explicit, so that you can
* write <code>Accumulator<double> a = 5;</code>.
*
* @param[in] y set \e sum = \e y.
**********************************************************************/
Accumulator(T y = T(0)) : _s(y), _t(0) {
static_assert(!std::numeric_limits<T>::is_integer,
"Accumulator type is not floating point");
}
/**
* Set the accumulator to a number.
*
* @param[in] y set \e sum = \e y.
**********************************************************************/
Accumulator& operator=(T y) { _s = y; _t = 0; return *this; }
/**
* Return the value held in the accumulator.
*
* @return \e sum.
**********************************************************************/
T operator()() const { return _s; }
/**
* Return the result of adding a number to \e sum (but don't change \e
* sum).
*
* @param[in] y the number to be added to the sum.
* @return \e sum + \e y.
**********************************************************************/
T operator()(T y) const { return Sum(y); }
/**
* Add a number to the accumulator.
*
* @param[in] y set \e sum += \e y.
**********************************************************************/
Accumulator& operator+=(T y) { Add(y); return *this; }
/**
* Subtract a number from the accumulator.
*
* @param[in] y set \e sum -= \e y.
**********************************************************************/
Accumulator& operator-=(T y) { Add(-y); return *this; }
/**
* Multiply accumulator by an integer. To avoid loss of accuracy, use only
* integers such that \e n &times; \e T is exactly representable as a \e T
* (i.e., &plusmn; powers of two). Use \e n = &minus;1 to negate \e sum.
*
* @param[in] n set \e sum *= \e n.
**********************************************************************/
Accumulator& operator*=(int n) { _s *= n; _t *= n; return *this; }
/**
* Multiply accumulator by a number. The fma (fused multiply and add)
* instruction is used (if available) in order to maintain accuracy.
*
* @param[in] y set \e sum *= \e y.
**********************************************************************/
Accumulator& operator*=(T y) {
using std::fma;
T d = _s; _s *= y;
d = fma(y, d, -_s); // the error in the first multiplication
_t = fma(y, _t, d); // add error to the second term
return *this;
}
/**
* Reduce accumulator to the range [-y/2, y/2].
*
* @param[in] y the modulus.
**********************************************************************/
Accumulator& remainder(T y) {
using std::remainder;
_s = remainder(_s, y);
Add(0); // This renormalizes the result.
return *this;
}
/**
* Test equality of an Accumulator with a number.
**********************************************************************/
bool operator==(T y) const { return _s == y; }
/**
* Test inequality of an Accumulator with a number.
**********************************************************************/
bool operator!=(T y) const { return _s != y; }
/**
* Less operator on an Accumulator and a number.
**********************************************************************/
bool operator<(T y) const { return _s < y; }
/**
* Less or equal operator on an Accumulator and a number.
**********************************************************************/
bool operator<=(T y) const { return _s <= y; }
/**
* Greater operator on an Accumulator and a number.
**********************************************************************/
bool operator>(T y) const { return _s > y; }
/**
* Greater or equal operator on an Accumulator and a number.
**********************************************************************/
bool operator>=(T y) const { return _s >= y; }
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_ACCUMULATOR_HPP

View File

@@ -0,0 +1,321 @@
/**
* \file AlbersEqualArea.hpp
* \brief Header for GeographicLib::AlbersEqualArea class
*
* Copyright (c) Charles Karney (2010-2021) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
#define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP 1
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief Albers equal area conic projection
*
* Implementation taken from the report,
* - J. P. Snyder,
* <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A
* Working Manual</a>, USGS Professional Paper 1395 (1987),
* pp. 101--102.
*
* This is a implementation of the equations in Snyder except that divided
* differences will be [have been] used to transform the expressions into
* ones which may be evaluated accurately. [In this implementation, the
* projection correctly becomes the cylindrical equal area or the azimuthal
* equal area projection when the standard latitude is the equator or a
* pole.]
*
* The ellipsoid parameters, the standard parallels, and the scale on the
* standard parallels are set in the constructor. Internally, the case with
* two standard parallels is converted into a single standard parallel, the
* latitude of minimum azimuthal scale, with an azimuthal scale specified on
* this parallel. This latitude is also used as the latitude of origin which
* is returned by AlbersEqualArea::OriginLatitude. The azimuthal scale on
* the latitude of origin is given by AlbersEqualArea::CentralScale. The
* case with two standard parallels at opposite poles is singular and is
* disallowed. The central meridian (which is a trivial shift of the
* longitude) is specified as the \e lon0 argument of the
* AlbersEqualArea::Forward and AlbersEqualArea::Reverse functions.
* AlbersEqualArea::Forward and AlbersEqualArea::Reverse also return the
* meridian convergence, &gamma;, and azimuthal scale, \e k. A small square
* aligned with the cardinal directions is projected to a rectangle with
* dimensions \e k (in the E-W direction) and 1/\e k (in the N-S direction).
* The E-W sides of the rectangle are oriented &gamma; degrees
* counter-clockwise from the \e x axis. There is no provision in this class
* for specifying a false easting or false northing or a different latitude
* of origin.
*
* Example of use:
* \include example-AlbersEqualArea.cpp
*
* <a href="ConicProj.1.html">ConicProj</a> is a command-line utility
* providing access to the functionality of LambertConformalConic and
* AlbersEqualArea.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT AlbersEqualArea {
private:
typedef Math::real real;
real eps_, epsx_, epsx2_, tol_, tol0_;
real _a, _f, _fm, _e2, _e, _e2m, _qZ, _qx;
real _sign, _lat0, _k0;
real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0;
static const int numit_ = 5; // Newton iterations in Reverse
static const int numit0_ = 20; // Newton iterations in Init
static real hyp(real x) {
using std::hypot;
return hypot(real(1), x);
}
// atanh( e * x)/ e if f > 0
// atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
// x if f = 0
real atanhee(real x) const {
using std::atan; using std::abs; using std::atanh;
return _f > 0 ? atanh(_e * x)/_e : (_f < 0 ? (atan(_e * x)/_e) : x);
}
// return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
static real atanhxm1(real x);
// Divided differences
// Definition: Df(x,y) = (f(x)-f(y))/(x-y)
// See:
// W. M. Kahan and R. J. Fateman,
// Symbolic computation of divided differences,
// SIGSAM Bull. 33(3), 7-28 (1999)
// https://doi.org/10.1145/334714.334716
// http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
//
// General rules
// h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
// h(x) = f(x)*g(x):
// Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y)
// = Df(x,y)*g(y) + Dg(x,y)*f(x)
// = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
//
// sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
static real Dsn(real x, real y, real sx, real sy) {
// sx = x/hyp(x)
real t = x * y;
return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) :
(x - y != 0 ? (sx - sy) / (x - y) : 1);
}
// Datanhee(x,y) = (atanee(x)-atanee(y))/(x-y)
// = atanhee((x-y)/(1-e^2*x*y))/(x-y)
real Datanhee(real x, real y) const {
real t = x - y, d = 1 - _e2 * x * y;
return t == 0 ? 1 / d :
(x*y < 0 ? atanhee(x) - atanhee(y) : atanhee(t / d)) / t;
}
// DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
real DDatanhee(real x, real y) const;
real DDatanhee0(real x, real y) const;
real DDatanhee1(real x, real y) const;
real DDatanhee2(real x, real y) const;
void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1);
real txif(real tphi) const;
real tphif(real txi) const;
friend class Ellipsoid; // For access to txif, tphif, etc.
public:
/**
* Constructor with a single standard parallel.
*
* @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @param[in] stdlat standard parallel (degrees), the circle of tangency.
* @param[in] k0 azimuthal scale on the standard parallel.
* @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k0 is
* not positive.
* @exception GeographicErr if \e stdlat is not in [&minus;90&deg;,
* 90&deg;].
**********************************************************************/
AlbersEqualArea(real a, real f, real stdlat, real k0);
/**
* Constructor with two standard parallels.
*
* @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @param[in] stdlat1 first standard parallel (degrees).
* @param[in] stdlat2 second standard parallel (degrees).
* @param[in] k1 azimuthal scale on the standard parallels.
* @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
* not positive.
* @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
* [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
* opposite poles.
**********************************************************************/
AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
/**
* Constructor with two standard parallels specified by sines and cosines.
*
* @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @param[in] sinlat1 sine of first standard parallel.
* @param[in] coslat1 cosine of first standard parallel.
* @param[in] sinlat2 sine of second standard parallel.
* @param[in] coslat2 cosine of second standard parallel.
* @param[in] k1 azimuthal scale on the standard parallels.
* @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k1 is
* not positive.
* @exception GeographicErr if \e stdlat1 or \e stdlat2 is not in
* [&minus;90&deg;, 90&deg;], or if \e stdlat1 and \e stdlat2 are
* opposite poles.
*
* This allows parallels close to the poles to be specified accurately.
* This routine computes the latitude of origin and the azimuthal scale at
* this latitude. If \e dlat = abs(\e lat2 &minus; \e lat1) &le; 160&deg;,
* then the error in the latitude of origin is less than 4.5 &times;
* 10<sup>&minus;14</sup>d;.
**********************************************************************/
AlbersEqualArea(real a, real f,
real sinlat1, real coslat1,
real sinlat2, real coslat2,
real k1);
/**
* Set the azimuthal scale for the projection.
*
* @param[in] lat (degrees).
* @param[in] k azimuthal scale at latitude \e lat (default 1).
* @exception GeographicErr \e k is not positive.
* @exception GeographicErr if \e lat is not in (&minus;90&deg;,
* 90&deg;).
*
* This allows a "latitude of conformality" to be specified.
**********************************************************************/
void SetScale(real lat, real k = real(1));
/**
* Forward projection, from geographic to Lambert conformal conic.
*
* @param[in] lon0 central meridian longitude (degrees).
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] gamma meridian convergence at point (degrees).
* @param[out] k azimuthal scale of projection at point; the radial
* scale is the 1/\e k.
*
* The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
* false easting or northing is added and \e lat should be in the range
* [&minus;90&deg;, 90&deg;]. The values of \e x and \e y returned for
* points which project to infinity (i.e., one or both of the poles) will
* be large but finite.
**********************************************************************/
void Forward(real lon0, real lat, real lon,
real& x, real& y, real& gamma, real& k) const;
/**
* Reverse projection, from Lambert conformal conic to geographic.
*
* @param[in] lon0 central meridian longitude (degrees).
* @param[in] x easting of point (meters).
* @param[in] y northing of point (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] gamma meridian convergence at point (degrees).
* @param[out] k azimuthal scale of projection at point; the radial
* scale is the 1/\e k.
*
* The latitude origin is given by AlbersEqualArea::LatitudeOrigin(). No
* false easting or northing is added. The value of \e lon returned is in
* the range [&minus;180&deg;, 180&deg;]. The value of \e lat returned is
* in the range [&minus;90&deg;, 90&deg;]. If the input point is outside
* the legal projected space the nearest pole is returned.
**********************************************************************/
void Reverse(real lon0, real x, real y,
real& lat, real& lon, real& gamma, real& k) const;
/**
* AlbersEqualArea::Forward without returning the convergence and
* scale.
**********************************************************************/
void Forward(real lon0, real lat, real lon,
real& x, real& y) const {
real gamma, k;
Forward(lon0, lat, lon, x, y, gamma, k);
}
/**
* AlbersEqualArea::Reverse without returning the convergence and
* scale.
**********************************************************************/
void Reverse(real lon0, real x, real y,
real& lat, real& lon) const {
real gamma, k;
Reverse(lon0, x, y, lat, lon, gamma, k);
}
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _a; }
/**
* @return \e f the flattening of the ellipsoid. This is the value used in
* the constructor.
**********************************************************************/
Math::real Flattening() const { return _f; }
/**
* @return latitude of the origin for the projection (degrees).
*
* This is the latitude of minimum azimuthal scale and equals the \e stdlat
* in the 1-parallel constructor and lies between \e stdlat1 and \e stdlat2
* in the 2-parallel constructors.
**********************************************************************/
Math::real OriginLatitude() const { return _lat0; }
/**
* @return central scale for the projection. This is the azimuthal scale
* on the latitude of origin.
**********************************************************************/
Math::real CentralScale() const { return _k0; }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
/**
* A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
* stdlat = 0, and \e k0 = 1. This degenerates to the cylindrical equal
* area projection.
**********************************************************************/
static const AlbersEqualArea& CylindricalEqualArea();
/**
* A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
* stdlat = 90&deg;, and \e k0 = 1. This degenerates to the
* Lambert azimuthal equal area projection.
**********************************************************************/
static const AlbersEqualArea& AzimuthalEqualAreaNorth();
/**
* A global instantiation of AlbersEqualArea with the WGS84 ellipsoid, \e
* stdlat = &minus;90&deg;, and \e k0 = 1. This degenerates to the
* Lambert azimuthal equal area projection.
**********************************************************************/
static const AlbersEqualArea& AzimuthalEqualAreaSouth();
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP

View File

@@ -0,0 +1,145 @@
/**
* \file AzimuthalEquidistant.hpp
* \brief Header for GeographicLib::AzimuthalEquidistant class
*
* Copyright (c) Charles Karney (2009-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP)
#define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP 1
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief Azimuthal equidistant projection
*
* Azimuthal equidistant projection centered at an arbitrary position on the
* ellipsoid. For a point in projected space (\e x, \e y), the geodesic
* distance from the center position is hypot(\e x, \e y) and the azimuth of
* the geodesic from the center point is atan2(\e x, \e y). The Forward and
* Reverse methods also return the azimuth \e azi of the geodesic at (\e x,
* \e y) and reciprocal scale \e rk in the azimuthal direction which,
* together with the basic properties of the projection, serve to specify
* completely the local affine transformation between geographic and
* projected coordinates.
*
* The conversions all take place using a Geodesic object (by default
* Geodesic::WGS84()). For more information on geodesics see \ref geodesic.
*
* Example of use:
* \include example-AzimuthalEquidistant.cpp
*
* <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility
* providing access to the functionality of AzimuthalEquidistant, Gnomonic,
* and CassiniSoldner.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT AzimuthalEquidistant {
private:
typedef Math::real real;
real eps_;
Geodesic _earth;
public:
/**
* Constructor for AzimuthalEquidistant.
*
* @param[in] earth the Geodesic object to use for geodesic calculations.
* By default this uses the WGS84 ellipsoid.
**********************************************************************/
explicit AzimuthalEquidistant(const Geodesic& earth = Geodesic::WGS84());
/**
* Forward projection, from geographic to azimuthal equidistant.
*
* @param[in] lat0 latitude of center point of projection (degrees).
* @param[in] lon0 longitude of center point of projection (degrees).
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] azi azimuth of geodesic at point (degrees).
* @param[out] rk reciprocal of azimuthal scale at point.
*
* \e lat0 and \e lat should be in the range [&minus;90&deg;, 90&deg;].
* The scale of the projection is 1 in the "radial" direction, \e azi
* clockwise from true north, and is 1/\e rk in the direction perpendicular
* to this. A call to Forward followed by a call to Reverse will return
* the original (\e lat, \e lon) (to within roundoff).
**********************************************************************/
void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y, real& azi, real& rk) const;
/**
* Reverse projection, from azimuthal equidistant to geographic.
*
* @param[in] lat0 latitude of center point of projection (degrees).
* @param[in] lon0 longitude of center point of projection (degrees).
* @param[in] x easting of point (meters).
* @param[in] y northing of point (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] azi azimuth of geodesic at point (degrees).
* @param[out] rk reciprocal of azimuthal scale at point.
*
* \e lat0 should be in the range [&minus;90&deg;, 90&deg;]. \e lat will
* be in the range [&minus;90&deg;, 90&deg;] and \e lon will be in the
* range [&minus;180&deg;, 180&deg;]. The scale of the projection is 1 in
* the "radial" direction, \e azi clockwise from true north, and is 1/\e rk
* in the direction perpendicular to this. A call to Reverse followed by a
* call to Forward will return the original (\e x, \e y) (to roundoff) only
* if the geodesic to (\e x, \e y) is a shortest path.
**********************************************************************/
void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon, real& azi, real& rk) const;
/**
* AzimuthalEquidistant::Forward without returning the azimuth and scale.
**********************************************************************/
void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y) const {
real azi, rk;
Forward(lat0, lon0, lat, lon, x, y, azi, rk);
}
/**
* AzimuthalEquidistant::Reverse without returning the azimuth and scale.
**********************************************************************/
void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon) const {
real azi, rk;
Reverse(lat0, lon0, x, y, lat, lon, azi, rk);
}
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }
/**
* @return \e f the flattening of the ellipsoid. This is the value
* inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real Flattening() const { return _earth.Flattening(); }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP

View File

@@ -0,0 +1,210 @@
/**
* \file CassiniSoldner.hpp
* \brief Header for GeographicLib::CassiniSoldner class
*
* Copyright (c) Charles Karney (2009-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP)
#define GEOGRAPHICLIB_CASSINISOLDNER_HPP 1
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp>
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief Cassini-Soldner projection
*
* Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e
* lon0, on the ellipsoid. This projection is a transverse cylindrical
* equidistant projection. The projection from (\e lat, \e lon) to easting
* and northing (\e x, \e y) is defined by geodesics as follows. Go north
* along a geodesic a distance \e y from the central point; then turn
* clockwise 90&deg; and go a distance \e x along a geodesic.
* (Although the initial heading is north, this changes to south if the pole
* is crossed.) This procedure uniquely defines the reverse projection. The
* forward projection is constructed as follows. Find the point (\e lat1, \e
* lon1) on the meridian closest to (\e lat, \e lon). Here we consider the
* full meridian so that \e lon1 may be either \e lon0 or \e lon0 +
* 180&deg;. \e x is the geodesic distance from (\e lat1, \e lon1) to
* (\e lat, \e lon), appropriately signed according to which side of the
* central meridian (\e lat, \e lon) lies. \e y is the shortest distance
* along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again,
* appropriately signed according to the initial heading. [Note that, in the
* case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e
* lon0) to (\e lat1, \e lon1) may not be the shortest path.] This procedure
* uniquely defines the forward projection except for a small class of points
* for which there may be two equally short routes for either leg of the
* path.
*
* Because of the properties of geodesics, the (\e x, \e y) grid is
* orthogonal. The scale in the easting direction is unity. The scale, \e
* k, in the northing direction is unity on the central meridian and
* increases away from the central meridian. The projection routines return
* \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the
* reciprocal of the scale in the northing direction.
*
* The conversions all take place using a Geodesic object (by default
* Geodesic::WGS84()). For more information on geodesics see \ref geodesic.
* The determination of (\e lat1, \e lon1) in the forward projection is by
* solving the inverse geodesic problem for (\e lat, \e lon) and its twin
* obtained by reflection in the meridional plane. The scale is found by
* determining where two neighboring geodesics intersecting the central
* meridian at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ratio
* of the reduced lengths for the two geodesics between that point and,
* respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
*
* Example of use:
* \include example-CassiniSoldner.cpp
*
* <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility
* providing access to the functionality of AzimuthalEquidistant, Gnomonic,
* and CassiniSoldner.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT CassiniSoldner {
private:
typedef Math::real real;
Geodesic _earth;
GeodesicLine _meridian;
real _sbet0, _cbet0;
static const unsigned maxit_ = 10;
public:
/**
* Constructor for CassiniSoldner.
*
* @param[in] earth the Geodesic object to use for geodesic calculations.
* By default this uses the WGS84 ellipsoid.
*
* This constructor makes an "uninitialized" object. Call Reset to set the
* central latitude and longitude, prior to calling Forward and Reverse.
**********************************************************************/
explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84());
/**
* Constructor for CassiniSoldner specifying a center point.
*
* @param[in] lat0 latitude of center point of projection (degrees).
* @param[in] lon0 longitude of center point of projection (degrees).
* @param[in] earth the Geodesic object to use for geodesic calculations.
* By default this uses the WGS84 ellipsoid.
*
* \e lat0 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
CassiniSoldner(real lat0, real lon0,
const Geodesic& earth = Geodesic::WGS84());
/**
* Set the central point of the projection
*
* @param[in] lat0 latitude of center point of projection (degrees).
* @param[in] lon0 longitude of center point of projection (degrees).
*
* \e lat0 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
void Reset(real lat0, real lon0);
/**
* Forward projection, from geographic to Cassini-Soldner.
*
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] azi azimuth of easting direction at point (degrees).
* @param[out] rk reciprocal of azimuthal northing scale at point.
*
* \e lat should be in the range [&minus;90&deg;, 90&deg;]. A call to
* Forward followed by a call to Reverse will return the original (\e lat,
* \e lon) (to within roundoff). The routine does nothing if the origin
* has not been set.
**********************************************************************/
void Forward(real lat, real lon,
real& x, real& y, real& azi, real& rk) const;
/**
* Reverse projection, from Cassini-Soldner to geographic.
*
* @param[in] x easting of point (meters).
* @param[in] y northing of point (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] azi azimuth of easting direction at point (degrees).
* @param[out] rk reciprocal of azimuthal northing scale at point.
*
* A call to Reverse followed by a call to Forward will return the original
* (\e x, \e y) (to within roundoff), provided that \e x and \e y are
* sufficiently small not to "wrap around" the earth. The routine does
* nothing if the origin has not been set.
**********************************************************************/
void Reverse(real x, real y,
real& lat, real& lon, real& azi, real& rk) const;
/**
* CassiniSoldner::Forward without returning the azimuth and scale.
**********************************************************************/
void Forward(real lat, real lon,
real& x, real& y) const {
real azi, rk;
Forward(lat, lon, x, y, azi, rk);
}
/**
* CassiniSoldner::Reverse without returning the azimuth and scale.
**********************************************************************/
void Reverse(real x, real y,
real& lat, real& lon) const {
real azi, rk;
Reverse(x, y, lat, lon, azi, rk);
}
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return true if the object has been initialized.
**********************************************************************/
bool Init() const { return _meridian.Init(); }
/**
* @return \e lat0 the latitude of origin (degrees).
**********************************************************************/
Math::real LatitudeOrigin() const
{ return _meridian.Latitude(); }
/**
* @return \e lon0 the longitude of origin (degrees).
**********************************************************************/
Math::real LongitudeOrigin() const
{ return _meridian.Longitude(); }
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }
/**
* @return \e f the flattening of the ellipsoid. This is the value
* inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real Flattening() const { return _earth.Flattening(); }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_CASSINISOLDNER_HPP

View File

@@ -0,0 +1,195 @@
/**
* \file CircularEngine.hpp
* \brief Header for GeographicLib::CircularEngine class
*
* Copyright (c) Charles Karney (2011-2015) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_CIRCULARENGINE_HPP)
#define GEOGRAPHICLIB_CIRCULARENGINE_HPP 1
#include <vector>
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp>
#if defined(_MSC_VER)
// Squelch warnings about dll vs vector
# pragma warning (push)
# pragma warning (disable: 4251)
#endif
namespace GeographicLib {
/**
* \brief Spherical harmonic sums for a circle
*
* The class is a companion to SphericalEngine. If the results of a
* spherical harmonic sum are needed for several points on a circle of
* constant latitude \e lat and height \e h, then SphericalEngine::Circle can
* compute the inner sum, which is independent of longitude \e lon, and
* produce a CircularEngine object. CircularEngine::operator()() can
* then be used to perform the outer sum for particular vales of \e lon.
* This can lead to substantial improvements in computational speed for high
* degree sum (approximately by a factor of \e N / 2 where \e N is the
* maximum degree).
*
* CircularEngine is tightly linked to the internals of SphericalEngine. For
* that reason, the constructor for this class is private. Use
* SphericalHarmonic::Circle, SphericalHarmonic1::Circle, and
* SphericalHarmonic2::Circle to create instances of this class.
*
* CircularEngine stores the coefficients needed to allow the summation over
* order to be performed in 2 or 6 vectors of length \e M + 1 (depending on
* whether gradients are to be calculated). For this reason the constructor
* may throw a std::bad_alloc exception.
*
* Example of use:
* \include example-CircularEngine.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT CircularEngine {
private:
typedef Math::real real;
enum normalization {
FULL = SphericalEngine::FULL,
SCHMIDT = SphericalEngine::SCHMIDT,
};
int _M;
bool _gradp;
unsigned _norm;
real _a, _r, _u, _t;
std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts;
real _q, _uq, _uq2;
Math::real Value(bool gradp, real sl, real cl,
real& gradx, real& grady, real& gradz) const;
friend class SphericalEngine;
CircularEngine(int M, bool gradp, unsigned norm,
real a, real r, real u, real t)
: _M(M)
, _gradp(gradp)
, _norm(norm)
, _a(a)
, _r(r)
, _u(u)
, _t(t)
, _wc(std::vector<real>(_M + 1, 0))
, _ws(std::vector<real>(_M + 1, 0))
, _wrc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wrs(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wtc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wts(std::vector<real>(_gradp ? _M + 1 : 0, 0))
{
_q = _a / _r;
_uq = _u * _q;
_uq2 = Math::sq(_uq);
}
void SetCoeff(int m, real wc, real ws)
{ _wc[m] = wc; _ws[m] = ws; }
void SetCoeff(int m, real wc, real ws,
real wrc, real wrs, real wtc, real wts) {
_wc[m] = wc; _ws[m] = ws;
if (_gradp) {
_wrc[m] = wrc; _wrs[m] = wrs;
_wtc[m] = wtc; _wts[m] = wts;
}
}
public:
/**
* A default constructor. CircularEngine::operator()() on the resulting
* object returns zero. The resulting object can be assigned to the result
* of SphericalHarmonic::Circle.
**********************************************************************/
CircularEngine()
: _M(-1)
, _gradp(true)
, _u(0)
, _t(1)
{}
/**
* Evaluate the sum for a particular longitude given in terms of its
* sine and cosine.
*
* @param[in] sinlon the sine of the longitude.
* @param[in] coslon the cosine of the longitude.
* @return \e V the value of the sum.
*
* The arguments must satisfy <i>sinlon</i><sup>2</sup> +
* <i>coslon</i><sup>2</sup> = 1.
**********************************************************************/
Math::real operator()(real sinlon, real coslon) const {
real dummy;
return Value(false, sinlon, coslon, dummy, dummy, dummy);
}
/**
* Evaluate the sum for a particular longitude.
*
* @param[in] lon the longitude (degrees).
* @return \e V the value of the sum.
**********************************************************************/
Math::real operator()(real lon) const {
real sinlon, coslon;
Math::sincosd(lon, sinlon, coslon);
return (*this)(sinlon, coslon);
}
/**
* Evaluate the sum and its gradient for a particular longitude given in
* terms of its sine and cosine.
*
* @param[in] sinlon the sine of the longitude.
* @param[in] coslon the cosine of the longitude.
* @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient.
* @return \e V the value of the sum.
*
* The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched. The arguments must satisfy <i>sinlon</i><sup>2</sup> +
* <i>coslon</i><sup>2</sup> = 1.
**********************************************************************/
Math::real operator()(real sinlon, real coslon,
real& gradx, real& grady, real& gradz) const {
return Value(true, sinlon, coslon, gradx, grady, gradz);
}
/**
* Evaluate the sum and its gradient for a particular longitude.
*
* @param[in] lon the longitude (degrees).
* @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient.
* @return \e V the value of the sum.
*
* The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched.
**********************************************************************/
Math::real operator()(real lon,
real& gradx, real& grady, real& gradz) const {
real sinlon, coslon;
Math::sincosd(lon, sinlon, coslon);
return (*this)(sinlon, coslon, gradx, grady, gradz);
}
};
} // namespace GeographicLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // GEOGRAPHICLIB_CIRCULARENGINE_HPP

25
external/include/GeographicLib/Config.h vendored Normal file
View File

@@ -0,0 +1,25 @@
#define GEOGRAPHICLIB_VERSION_STRING "1.52"
#define GEOGRAPHICLIB_VERSION_MAJOR 1
#define GEOGRAPHICLIB_VERSION_MINOR 52
#define GEOGRAPHICLIB_VERSION_PATCH 0
#define GEOGRAPHICLIB_DATA "C:/Users/sven/Documents/Visual Studio 2019/Projects/TST/install"
// These are macros which affect the building of the library
#define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
#define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
#define GEOGRAPHICLIB_PRECISION 1
// Specify whether GeographicLib is a shared or static library. When compiling
// under Visual Studio it is necessary to specify whether GeographicLib is a
// shared library. This is done with the macro GEOGRAPHICLIB_SHARED_LIB, which
// cmake will correctly define as 0 or 1 when only one type of library is in
// the package. If both shared and static libraries are available,
// GEOGRAPHICLIB_SHARED_LIB is set to 2 which triggers a preprocessor error in
// Constants.hpp. In this case, the appropriate value (0 or 1) for
// GEOGRAPHICLIB_SHARED_LIB must be specified when compiling any program that
// includes GeographicLib headers. This is done automatically if GeographicLib
// and the user's code were built with cmake version 2.8.11 (which introduced
// the command target_compile_definitions) or later.
#if !defined(GEOGRAPHICLIB_SHARED_LIB)
#define GEOGRAPHICLIB_SHARED_LIB 0
#endif

View File

@@ -0,0 +1,329 @@
/**
* \file Constants.hpp
* \brief Header for GeographicLib::Constants class
*
* Copyright (c) Charles Karney (2008-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_CONSTANTS_HPP)
#define GEOGRAPHICLIB_CONSTANTS_HPP 1
#include <GeographicLib/Config.h>
/**
* @relates GeographicLib::Constants
* Pack the version components into a single integer. Users should not rely on
* this particular packing of the components of the version number; see the
* documentation for GEOGRAPHICLIB_VERSION, below.
**********************************************************************/
#define GEOGRAPHICLIB_VERSION_NUM(a,b,c) ((((a) * 10000 + (b)) * 100) + (c))
/**
* @relates GeographicLib::Constants
* The version of GeographicLib as a single integer, packed as MMmmmmpp where
* MM is the major version, mmmm is the minor version, and pp is the patch
* level. Users should not rely on this particular packing of the components
* of the version number. Instead they should use a test such as \code
#if GEOGRAPHICLIB_VERSION >= GEOGRAPHICLIB_VERSION_NUM(1,37,0)
...
#endif
* \endcode
**********************************************************************/
#define GEOGRAPHICLIB_VERSION \
GEOGRAPHICLIB_VERSION_NUM(GEOGRAPHICLIB_VERSION_MAJOR, \
GEOGRAPHICLIB_VERSION_MINOR, \
GEOGRAPHICLIB_VERSION_PATCH)
// For reference, here is a table of Visual Studio and _MSC_VER
// correspondences:
//
// _MSC_VER Visual Studio
// 1100 vc5
// 1200 vc6
// 1300 vc7
// 1310 vc7.1 (2003)
// 1400 vc8 (2005)
// 1500 vc9 (2008)
// 1600 vc10 (2010)
// 1700 vc11 (2012)
// 1800 vc12 (2013)
// 1900 vc14 (2015) First version of VS to include enough C++11 support
// 191[0-9] vc15 (2017)
// 192[0-9] vc16 (2019)
#if defined(_MSC_VER) && defined(GEOGRAPHICLIB_SHARED_LIB) && \
GEOGRAPHICLIB_SHARED_LIB
# if GEOGRAPHICLIB_SHARED_LIB > 1
# error GEOGRAPHICLIB_SHARED_LIB must be 0 or 1
# elif defined(GeographicLib_SHARED_EXPORTS)
# define GEOGRAPHICLIB_EXPORT __declspec(dllexport)
# else
# define GEOGRAPHICLIB_EXPORT __declspec(dllimport)
# endif
#else
# define GEOGRAPHICLIB_EXPORT
#endif
// Use GEOGRAPHICLIB_DEPRECATED to mark functions, types or variables as
// deprecated. Code inspired by Apache Subversion's svn_types.h file (via
// MPFR).
#if defined(__GNUC__)
# if __GNUC__ > 4
# define GEOGRAPHICLIB_DEPRECATED(msg) __attribute__((deprecated(msg)))
# else
# define GEOGRAPHICLIB_DEPRECATED(msg) __attribute__((deprecated))
# endif
#elif defined(_MSC_VER) && _MSC_VER >= 1300
# define GEOGRAPHICLIB_DEPRECATED(msg) __declspec(deprecated(msg))
#else
# define GEOGRAPHICLIB_DEPRECATED(msg)
#endif
#include <stdexcept>
#include <string>
#include <GeographicLib/Math.hpp>
/**
* \brief Namespace for %GeographicLib
*
* All of %GeographicLib is defined within the GeographicLib namespace. In
* addition all the header files are included via %GeographicLib/Class.hpp.
* This minimizes the likelihood of conflicts with other packages.
**********************************************************************/
namespace GeographicLib {
/**
* \brief %Constants needed by %GeographicLib
*
* Define constants specifying the WGS84 ellipsoid, the UTM and UPS
* projections, and various unit conversions.
*
* Example of use:
* \include example-Constants.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT Constants {
private:
typedef Math::real real;
Constants(); // Disable constructor
public:
/**
* A synonym for Math::degree<real>().
**********************************************************************/
static Math::real degree() { return Math::degree(); }
/**
* @return the number of radians in an arcminute.
**********************************************************************/
static Math::real arcminute()
{ return Math::degree() / 60; }
/**
* @return the number of radians in an arcsecond.
**********************************************************************/
static Math::real arcsecond()
{ return Math::degree() / 3600; }
/** \name Ellipsoid parameters
**********************************************************************/
///@{
/**
* @tparam T the type of the returned value.
* @return the equatorial radius of WGS84 ellipsoid (6378137 m).
**********************************************************************/
template<typename T = real> static T WGS84_a()
{ return 6378137 * meter<T>(); }
/**
* @tparam T the type of the returned value.
* @return the flattening of WGS84 ellipsoid (1/298.257223563).
**********************************************************************/
template<typename T = real> static T WGS84_f() {
// Evaluating this as 1000000000 / T(298257223563LL) reduces the
// round-off error by about 10%. However, expressing the flattening as
// 1/298.257223563 is well ingrained.
return 1 / ( T(298257223563LL) / 1000000000 );
}
/**
* @tparam T the type of the returned value.
* @return the gravitational constant of the WGS84 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>&minus;2</sup>.
**********************************************************************/
template<typename T = real> static T WGS84_GM()
{ return T(3986004) * 100000000 + 41800000; }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the WGS84 ellipsoid, &omega;, in rad
* s<sup>&minus;1</sup>.
**********************************************************************/
template<typename T = real> static T WGS84_omega()
{ return 7292115 / (T(1000000) * 100000); }
/**
* @tparam T the type of the returned value.
* @return the equatorial radius of GRS80 ellipsoid, \e a, in m.
**********************************************************************/
template<typename T = real> static T GRS80_a()
{ return 6378137 * meter<T>(); }
/**
* @tparam T the type of the returned value.
* @return the gravitational constant of the GRS80 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>&minus;2</sup>.
**********************************************************************/
template<typename T = real> static T GRS80_GM()
{ return T(3986005) * 100000000; }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the GRS80 ellipsoid, &omega;, in rad
* s<sup>&minus;1</sup>.
*
* This is about 2 &pi; 366.25 / (365.25 &times; 24 &times; 3600) rad
* s<sup>&minus;1</sup>. 365.25 is the number of days in a Julian year and
* 365.35/366.25 converts from solar days to sidereal days. Using the
* number of days in a Gregorian year (365.2425) results in a worse
* approximation (because the Gregorian year includes the precession of the
* earth's axis).
**********************************************************************/
template<typename T = real> static T GRS80_omega()
{ return 7292115 / (T(1000000) * 100000); }
/**
* @tparam T the type of the returned value.
* @return the dynamical form factor of the GRS80 ellipsoid,
* <i>J</i><sub>2</sub>.
**********************************************************************/
template<typename T = real> static T GRS80_J2()
{ return T(108263) / 100000000; }
/**
* @tparam T the type of the returned value.
* @return the central scale factor for UTM (0.9996).
**********************************************************************/
template<typename T = real> static T UTM_k0()
{return T(9996) / 10000; }
/**
* @tparam T the type of the returned value.
* @return the central scale factor for UPS (0.994).
**********************************************************************/
template<typename T = real> static T UPS_k0()
{ return T(994) / 1000; }
///@}
/** \name SI units
**********************************************************************/
///@{
/**
* @tparam T the type of the returned value.
* @return the number of meters in a meter.
*
* This is unity, but this lets the internal system of units be changed if
* necessary.
**********************************************************************/
template<typename T = real> static T meter() { return T(1); }
/**
* @return the number of meters in a kilometer.
**********************************************************************/
static Math::real kilometer()
{ return 1000 * meter<real>(); }
/**
* @return the number of meters in a nautical mile (approximately 1 arc
* minute)
**********************************************************************/
static Math::real nauticalmile()
{ return 1852 * meter<real>(); }
/**
* @tparam T the type of the returned value.
* @return the number of square meters in a square meter.
*
* This is unity, but this lets the internal system of units be changed if
* necessary.
**********************************************************************/
template<typename T = real> static T square_meter()
{ return meter<T>() * meter<T>(); }
/**
* @return the number of square meters in a hectare.
**********************************************************************/
static Math::real hectare()
{ return 10000 * square_meter<real>(); }
/**
* @return the number of square meters in a square kilometer.
**********************************************************************/
static Math::real square_kilometer()
{ return kilometer() * kilometer(); }
/**
* @return the number of square meters in a square nautical mile.
**********************************************************************/
static Math::real square_nauticalmile()
{ return nauticalmile() * nauticalmile(); }
///@}
/** \name Anachronistic British units
**********************************************************************/
///@{
/**
* @return the number of meters in an international foot.
**********************************************************************/
static Math::real foot()
{ return real(254 * 12) / 10000 * meter<real>(); }
/**
* @return the number of meters in a yard.
**********************************************************************/
static Math::real yard() { return 3 * foot(); }
/**
* @return the number of meters in a fathom.
**********************************************************************/
static Math::real fathom() { return 2 * yard(); }
/**
* @return the number of meters in a chain.
**********************************************************************/
static Math::real chain() { return 22 * yard(); }
/**
* @return the number of meters in a furlong.
**********************************************************************/
static Math::real furlong() { return 10 * chain(); }
/**
* @return the number of meters in a statute mile.
**********************************************************************/
static Math::real mile() { return 8 * furlong(); }
/**
* @return the number of square meters in an acre.
**********************************************************************/
static Math::real acre() { return chain() * furlong(); }
/**
* @return the number of square meters in a square statute mile.
**********************************************************************/
static Math::real square_mile() { return mile() * mile(); }
///@}
/** \name Anachronistic US units
**********************************************************************/
///@{
/**
* @return the number of meters in a US survey foot.
**********************************************************************/
static Math::real surveyfoot()
{ return real(1200) / 3937 * meter<real>(); }
///@}
};
/**
* \brief Exception handling for %GeographicLib
*
* A class to handle exceptions. It's derived from std::runtime_error so it
* can be caught by the usual catch clauses.
*
* Example of use:
* \include example-GeographicErr.cpp
**********************************************************************/
class GeographicErr : public std::runtime_error {
public:
/**
* Constructor
*
* @param[in] msg a string message, which is accessible in the catch
* clause via what().
**********************************************************************/
GeographicErr(const std::string& msg) : std::runtime_error(msg) {}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_CONSTANTS_HPP

405
external/include/GeographicLib/DMS.hpp vendored Normal file
View File

@@ -0,0 +1,405 @@
/**
* \file DMS.hpp
* \brief Header for GeographicLib::DMS class
*
* Copyright (c) Charles Karney (2008-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_DMS_HPP)
#define GEOGRAPHICLIB_DMS_HPP 1
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Utility.hpp>
#if defined(_MSC_VER)
// Squelch warnings about dll vs vector and constant conditional expressions
# pragma warning (push)
# pragma warning (disable: 4251 4127)
#endif
namespace GeographicLib {
/**
* \brief Convert between degrees and the %DMS representation
*
* Parse a string representing degree, minutes, and seconds and return the
* angle in degrees and format an angle in degrees as degree, minutes, and
* seconds. In addition, handle NANs and infinities on input and output.
*
* Example of use:
* \include example-DMS.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT DMS {
public:
/**
* Indicator for presence of hemisphere indicator (N/S/E/W) on latitudes
* and longitudes.
**********************************************************************/
enum flag {
/**
* No indicator present.
* @hideinitializer
**********************************************************************/
NONE = 0,
/**
* Latitude indicator (N/S) present.
* @hideinitializer
**********************************************************************/
LATITUDE = 1,
/**
* Longitude indicator (E/W) present.
* @hideinitializer
**********************************************************************/
LONGITUDE = 2,
/**
* Used in Encode to indicate output of an azimuth in [000, 360) with no
* letter indicator.
* @hideinitializer
**********************************************************************/
AZIMUTH = 3,
/**
* Used in Encode to indicate output of a plain number.
* @hideinitializer
**********************************************************************/
NUMBER = 4,
};
/**
* Indicator for trailing units on an angle.
**********************************************************************/
enum component {
/**
* Trailing unit is degrees.
* @hideinitializer
**********************************************************************/
DEGREE = 0,
/**
* Trailing unit is arc minutes.
* @hideinitializer
**********************************************************************/
MINUTE = 1,
/**
* Trailing unit is arc seconds.
* @hideinitializer
**********************************************************************/
SECOND = 2,
};
private:
typedef Math::real real;
// Replace all occurrences of pat by c. If c is NULL remove pat.
static void replace(std::string& s, const std::string& pat, char c) {
std::string::size_type p = 0;
int count = c ? 1 : 0;
while (true) {
p = s.find(pat, p);
if (p == std::string::npos)
break;
s.replace(p, pat.length(), count, c);
}
}
static const char* const hemispheres_;
static const char* const signs_;
static const char* const digits_;
static const char* const dmsindicators_;
static const char* const components_[3];
static Math::real NumMatch(const std::string& s);
static Math::real InternalDecode(const std::string& dmsa, flag& ind);
DMS(); // Disable constructor
public:
/**
* Convert a string in DMS to an angle.
*
* @param[in] dms string input.
* @param[out] ind a DMS::flag value signaling the presence of a
* hemisphere indicator.
* @exception GeographicErr if \e dms is malformed (see below).
* @return angle (degrees).
*
* Degrees, minutes, and seconds are indicated by the characters d, '
* (single quote), &quot; (double quote), and these components may only be
* given in this order. Any (but not all) components may be omitted and
* other symbols (e.g., the &deg; symbol for degrees and the unicode prime
* and double prime symbols for minutes and seconds) may be substituted;
* two single quotes can be used instead of &quot;. The last component
* indicator may be omitted and is assumed to be the next smallest unit
* (thus 33d10 is interpreted as 33d10'). The final component may be a
* decimal fraction but the non-final components must be integers. Instead
* of using d, ', and &quot; to indicate degrees, minutes, and seconds, :
* (colon) may be used to <i>separate</i> these components (numbers must
* appear before and after each colon); thus 50d30'10.3&quot; may be
* written as 50:30:10.3, 5.5' may be written 0:5.5, and so on. The
* integer parts of the minutes and seconds components must be less
* than 60. A single leading sign is permitted. A hemisphere designator
* (N, E, W, S) may be added to the beginning or end of the string. The
* result is multiplied by the implied sign of the hemisphere designator
* (negative for S and W). In addition \e ind is set to DMS::LATITUDE if N
* or S is present, to DMS::LONGITUDE if E or W is present, and to
* DMS::NONE otherwise. Throws an error on a malformed string. No check
* is performed on the range of the result. Examples of legal and illegal
* strings are
* - <i>LEGAL</i> (all the entries on each line are equivalent)
* - -20.51125, 20d30'40.5&quot;S, -20&deg;30'40.5, -20d30.675,
* N-20d30'40.5&quot;, -20:30:40.5
* - 4d0'9, 4d9&quot;, 4d9'', 4:0:9, 004:00:09, 4.0025, 4.0025d, 4d0.15,
* 04:.15
* - 4:59.99999999999999, 4:60.0, 4:59:59.9999999999999, 4:59:60.0, 5
* - <i>ILLEGAL</i> (the exception thrown explains the problem)
* - 4d5&quot;4', 4::5, 4:5:, :4:5, 4d4.5'4&quot;, -N20.5, 1.8e2d, 4:60,
* 4:59:60
*
* The decoding operation can also perform addition and subtraction
* operations. If the string includes <i>internal</i> signs (i.e., not at
* the beginning nor immediately after an initial hemisphere designator),
* then the string is split immediately before such signs and each piece is
* decoded according to the above rules and the results added; thus
* <code>S3-2.5+4.1N</code> is parsed as the sum of <code>S3</code>,
* <code>-2.5</code>, <code>+4.1N</code>. Any piece can include a
* hemisphere designator; however, if multiple designators are given, they
* must compatible; e.g., you cannot mix N and E. In addition, the
* designator can appear at the beginning or end of the first piece, but
* must be at the end of all subsequent pieces (a hemisphere designator is
* not allowed after the initial sign). Examples of legal and illegal
* combinations are
* - <i>LEGAL</i> (these are all equivalent)
* - 070:00:45, 70:01:15W+0:0.5, 70:01:15W-0:0:30W, W70:01:15+0:0:30E
* - <i>ILLEGAL</i> (the exception thrown explains the problem)
* - 70:01:15W+0:0:15N, W70:01:15+W0:0:15
*
* \warning The "exponential" notation is not recognized. Thus
* <code>7.0E1</code> is illegal, while <code>7.0E+1</code> is parsed as
* <code>(7.0E) + (+1)</code>, yielding the same result as
* <code>8.0E</code>.
*
* \note At present, all the string handling in the C++ implementation of
* %GeographicLib is with 8-bit characters. The support for unicode
* symbols for degrees, minutes, and seconds is therefore via the
* <a href="https://en.wikipedia.org/wiki/UTF-8">UTF-8</a> encoding. (The
* JavaScript implementation of this class uses unicode natively, of
* course.)
*
* Here is the list of Unicode symbols supported for degrees, minutes,
* seconds, and the plus and minus signs; various symbols denoting variants
* of a space, which may separate the components of a DMS string, are
* removed:
* - degrees:
* - d, D lower and upper case letters
* - U+00b0 degree symbol (&deg;)
* - U+00ba masculine ordinal indicator (&ordm;)
* - U+2070 superscript zero (⁰)
* - U+02da ring above (˚)
* - U+2218 compose function (∘)
* - * the <a href="https://grid.nga.mil">GRiD</a> symbol for degrees
* - minutes:
* - ' apostrophe
* - ` grave accent
* - U+2032 prime (&prime;)
* - U+2035 back prime ()
* - U+00b4 acute accent (&acute;)
* - U+2018 left single quote (&lsquo;)
* - U+2019 right single quote (&rsquo;)
* - U+201b reversed-9 single quote ()
* - U+02b9 modifier letter prime (ʹ)
* - U+02ca modifier letter acute accent (ˊ)
* - U+02cb modifier letter grave accent (ˋ)
* - seconds:
* - &quot; quotation mark
* - U+2033 double prime (&Prime;)
* - U+2036 reversed double prime (‶)
* + U+02dd double acute accent (˝)
* - U+201c left double quote (&ldquo;)
* - U+201d right double quote (&rdquo;)
* - U+201f reversed-9 double quote (‟)
* - U+02ba modifier letter double prime (ʺ)
* - '&nbsp;' any two consecutive symbols for minutes
* - plus sign:
* - + plus
* - U+2795 heavy plus ()
* - U+2064 invisible plus (||)
* - minus sign:
* - - hyphen
* - U+2010 dash ()
* - U+2011 non-breaking hyphen ()
* - U+2013 en dash (&ndash;)
* - U+2014 em dash (&mdash;)
* - U+2212 minus sign (&minus;)
* - U+2796 heavy minus ()
* - ignored spaces:
* - U+00a0 non-breaking space
* - U+2007 figure space (||)
* - U+2009 thin space (|&thinsp;|)
* - U+200a hair space (||)
* - U+200b invisible space (||)
* - U+202f narrow space (||)
* - U+2063 invisible separator (||)
* .
* The codes with a leading zero byte, e.g., U+00b0, are accepted in their
* UTF-8 coded form 0xc2 0xb0 and as a single byte 0xb0.
**********************************************************************/
static Math::real Decode(const std::string& dms, flag& ind);
/**
* Convert DMS to an angle.
*
* @param[in] d degrees.
* @param[in] m arc minutes.
* @param[in] s arc seconds.
* @return angle (degrees)
*
* This does not propagate the sign on \e d to the other components,
* so -3d20' would need to be represented as - DMS::Decode(3.0, 20.0) or
* DMS::Decode(-3.0, -20.0).
**********************************************************************/
static Math::real Decode(real d, real m = 0, real s = 0)
{ return d + (m + s / 60) / 60; }
/**
* Convert a pair of strings to latitude and longitude.
*
* @param[in] dmsa first string.
* @param[in] dmsb second string.
* @param[out] lat latitude (degrees).
* @param[out] lon longitude (degrees).
* @param[in] longfirst if true assume longitude is given before latitude
* in the absence of hemisphere designators (default false).
* @exception GeographicErr if \e dmsa or \e dmsb is malformed.
* @exception GeographicErr if \e dmsa and \e dmsb are both interpreted as
* latitudes.
* @exception GeographicErr if \e dmsa and \e dmsb are both interpreted as
* longitudes.
* @exception GeographicErr if decoded latitude is not in [&minus;90&deg;,
* 90&deg;].
*
* By default, the \e lat (resp., \e lon) is assigned to the results of
* decoding \e dmsa (resp., \e dmsb). However this is overridden if either
* \e dmsa or \e dmsb contain a latitude or longitude hemisphere designator
* (N, S, E, W). If an exception is thrown, \e lat and \e lon are
* unchanged.
**********************************************************************/
static void DecodeLatLon(const std::string& dmsa, const std::string& dmsb,
real& lat, real& lon,
bool longfirst = false);
/**
* Convert a string to an angle in degrees.
*
* @param[in] angstr input string.
* @exception GeographicErr if \e angstr is malformed.
* @exception GeographicErr if \e angstr includes a hemisphere designator.
* @return angle (degrees)
*
* No hemisphere designator is allowed and no check is done on the range of
* the result.
**********************************************************************/
static Math::real DecodeAngle(const std::string& angstr);
/**
* Convert a string to an azimuth in degrees.
*
* @param[in] azistr input string.
* @exception GeographicErr if \e azistr is malformed.
* @exception GeographicErr if \e azistr includes a N/S designator.
* @return azimuth (degrees) reduced to the range [&minus;180&deg;,
* 180&deg;].
*
* A hemisphere designator E/W can be used; the result is multiplied by
* &minus;1 if W is present.
**********************************************************************/
static Math::real DecodeAzimuth(const std::string& azistr);
/**
* Convert angle (in degrees) into a DMS string (using d, ', and &quot;).
*
* @param[in] angle input angle (degrees)
* @param[in] trailing DMS::component value indicating the trailing units
* of the string (this component is given as a decimal number if
* necessary).
* @param[in] prec the number of digits after the decimal point for the
* trailing component.
* @param[in] ind DMS::flag value indicating additional formatting.
* @param[in] dmssep if non-null, use as the DMS separator character
* (instead of d, ', &quot; delimiters).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return formatted string
*
* The interpretation of \e ind is as follows:
* - ind == DMS::NONE, signed result no leading zeros on degrees except in
* the units place, e.g., -8d03'.
* - ind == DMS::LATITUDE, trailing N or S hemisphere designator, no sign,
* pad degrees to 2 digits, e.g., 08d03'S.
* - ind == DMS::LONGITUDE, trailing E or W hemisphere designator, no
* sign, pad degrees to 3 digits, e.g., 008d03'W.
* - ind == DMS::AZIMUTH, convert to the range [0, 360&deg;), no
* sign, pad degrees to 3 digits, e.g., 351d57'.
* .
* The integer parts of the minutes and seconds components are always given
* with 2 digits.
**********************************************************************/
static std::string Encode(real angle, component trailing, unsigned prec,
flag ind = NONE, char dmssep = char(0));
/**
* Convert angle into a DMS string (using d, ', and &quot;) selecting the
* trailing component based on the precision.
*
* @param[in] angle input angle (degrees)
* @param[in] prec the precision relative to 1 degree.
* @param[in] ind DMS::flag value indicated additional formatting.
* @param[in] dmssep if non-null, use as the DMS separator character
* (instead of d, ', &quot; delimiters).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return formatted string
*
* \e prec indicates the precision relative to 1 degree, e.g., \e prec = 3
* gives a result accurate to 0.1' and \e prec = 4 gives a result accurate
* to 1&quot;. \e ind is interpreted as in DMS::Encode with the additional
* facility that DMS::NUMBER represents \e angle as a number in fixed
* format with precision \e prec.
**********************************************************************/
static std::string Encode(real angle, unsigned prec, flag ind = NONE,
char dmssep = char(0)) {
return ind == NUMBER ? Utility::str(angle, int(prec)) :
Encode(angle,
prec < 2 ? DEGREE : (prec < 4 ? MINUTE : SECOND),
prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4),
ind, dmssep);
}
/**
* Split angle into degrees and minutes
*
* @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes.
**********************************************************************/
static void Encode(real ang, real& d, real& m) {
d = int(ang); m = 60 * (ang - d);
}
/**
* Split angle into degrees and minutes and seconds.
*
* @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes (an integer returned as a real)
* @param[out] s arc seconds.
**********************************************************************/
static void Encode(real ang, real& d, real& m, real& s) {
d = int(ang); ang = 60 * (ang - d);
m = int(ang); s = 60 * (ang - m);
}
};
} // namespace GeographicLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // GEOGRAPHICLIB_DMS_HPP

View File

@@ -0,0 +1,542 @@
/**
* \file Ellipsoid.hpp
* \brief Header for GeographicLib::Ellipsoid class
*
* Copyright (c) Charles Karney (2012-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_ELLIPSOID_HPP)
#define GEOGRAPHICLIB_ELLIPSOID_HPP 1
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/TransverseMercator.hpp>
#include <GeographicLib/EllipticFunction.hpp>
#include <GeographicLib/AlbersEqualArea.hpp>
namespace GeographicLib {
/**
* \brief Properties of an ellipsoid
*
* This class returns various properties of the ellipsoid and converts
* between various types of latitudes. The latitude conversions are also
* possible using the various projections supported by %GeographicLib; but
* Ellipsoid provides more direct access (sometimes using private functions
* of the projection classes). Ellipsoid::RectifyingLatitude,
* Ellipsoid::InverseRectifyingLatitude, and Ellipsoid::MeridianDistance
* provide functionality which can be provided by the Geodesic class.
* However Geodesic uses a series approximation (valid for abs \e f < 1/150),
* whereas Ellipsoid computes these quantities using EllipticFunction which
* provides accurate results even when \e f is large. Use of this class
* should be limited to &minus;3 < \e f < 3/4 (i.e., 1/4 < b/a < 4).
*
* Example of use:
* \include example-Ellipsoid.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT Ellipsoid {
private:
typedef Math::real real;
static const int numit_ = 10;
real stol_;
real _a, _f, _f1, _f12, _e2, _es, _e12, _n, _b;
TransverseMercator _tm;
EllipticFunction _ell;
AlbersEqualArea _au;
// These are the alpha and beta coefficients in the Krueger series from
// TransverseMercator. Thy are used by RhumbSolve to compute
// (psi2-psi1)/(mu2-mu1).
const Math::real* ConformalToRectifyingCoeffs() const { return _tm._alp; }
const Math::real* RectifyingToConformalCoeffs() const { return _tm._bet; }
friend class Rhumb; friend class RhumbLine;
public:
/** \name Constructor
**********************************************************************/
///@{
/**
* Constructor for a ellipsoid with
*
* @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
* positive.
**********************************************************************/
Ellipsoid(real a, real f);
///@}
/** \name %Ellipsoid dimensions.
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _a; }
/**
* @return \e b the polar semi-axis (meters).
**********************************************************************/
Math::real MinorRadius() const { return _b; }
/**
* @return \e L the distance between the equator and a pole along a
* meridian (meters). For a sphere \e L = (&pi;/2) \e a. The radius
* of a sphere with the same meridian length is \e L / (&pi;/2).
**********************************************************************/
Math::real QuarterMeridian() const;
/**
* @return \e A the total area of the ellipsoid (meters<sup>2</sup>). For
* a sphere \e A = 4&pi; <i>a</i><sup>2</sup>. The radius of a sphere
* with the same area is sqrt(\e A / (4&pi;)).
**********************************************************************/
Math::real Area() const;
/**
* @return \e V the total volume of the ellipsoid (meters<sup>3</sup>).
* For a sphere \e V = (4&pi; / 3) <i>a</i><sup>3</sup>. The radius of
* a sphere with the same volume is cbrt(\e V / (4&pi;/3)).
**********************************************************************/
Math::real Volume() const
{ return (4 * Math::pi()) * Math::sq(_a) * _b / 3; }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
/** \name %Ellipsoid shape
**********************************************************************/
///@{
/**
* @return \e f = (\e a &minus; \e b) / \e a, the flattening of the
* ellipsoid. This is the value used in the constructor. This is zero,
* positive, or negative for a sphere, oblate ellipsoid, or prolate
* ellipsoid.
**********************************************************************/
Math::real Flattening() const { return _f; }
/**
* @return \e f ' = (\e a &minus; \e b) / \e b, the second flattening of
* the ellipsoid. This is zero, positive, or negative for a sphere,
* oblate ellipsoid, or prolate ellipsoid.
**********************************************************************/
Math::real SecondFlattening() const { return _f / (1 - _f); }
/**
* @return \e n = (\e a &minus; \e b) / (\e a + \e b), the third flattening
* of the ellipsoid. This is zero, positive, or negative for a sphere,
* oblate ellipsoid, or prolate ellipsoid.
**********************************************************************/
Math::real ThirdFlattening() const { return _n; }
/**
* @return <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity squared
* of the ellipsoid. This is zero, positive, or negative for a sphere,
* oblate ellipsoid, or prolate ellipsoid.
**********************************************************************/
Math::real EccentricitySq() const { return _e2; }
/**
* @return <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity
* squared of the ellipsoid. This is zero, positive, or negative for a
* sphere, oblate ellipsoid, or prolate ellipsoid.
**********************************************************************/
Math::real SecondEccentricitySq() const { return _e12; }
/**
* @return <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> + <i>b</i><sup>2</sup>),
* the third eccentricity squared of the ellipsoid. This is zero,
* positive, or negative for a sphere, oblate ellipsoid, or prolate
* ellipsoid.
**********************************************************************/
Math::real ThirdEccentricitySq() const { return _e2 / (2 - _e2); }
///@}
/** \name Latitude conversion.
**********************************************************************/
///@{
/**
* @param[in] phi the geographic latitude (degrees).
* @return &beta; the parametric latitude (degrees).
*
* The geographic latitude, &phi;, is the angle between the equatorial
* plane and a vector normal to the surface of the ellipsoid.
*
* The parametric latitude (also called the reduced latitude), &beta;,
* allows the cartesian coordinated of a meridian to be expressed
* conveniently in parametric form as
* - \e R = \e a cos &beta;
* - \e Z = \e b sin &beta;
* .
* where \e a and \e b are the equatorial radius and the polar semi-axis.
* For a sphere &beta; = &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &beta; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real ParametricLatitude(real phi) const;
/**
* @param[in] beta the parametric latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* &beta; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &phi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real InverseParametricLatitude(real beta) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &theta; the geocentric latitude (degrees).
*
* The geocentric latitude, &theta;, is the angle between the equatorial
* plane and a line between the center of the ellipsoid and a point on the
* ellipsoid. For a sphere &theta; = &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &theta; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real GeocentricLatitude(real phi) const;
/**
* @param[in] theta the geocentric latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* &theta; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &phi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real InverseGeocentricLatitude(real theta) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &mu; the rectifying latitude (degrees).
*
* The rectifying latitude, &mu;, has the property that the distance along
* a meridian of the ellipsoid between two points with rectifying latitudes
* &mu;<sub>1</sub> and &mu;<sub>2</sub> is equal to
* (&mu;<sub>2</sub> - &mu;<sub>1</sub>) \e L / 90&deg;,
* where \e L = QuarterMeridian(). For a sphere &mu; = &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &mu; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real RectifyingLatitude(real phi) const;
/**
* @param[in] mu the rectifying latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* &mu; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &phi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real InverseRectifyingLatitude(real mu) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &xi; the authalic latitude (degrees).
*
* The authalic latitude, &xi;, has the property that the area of the
* ellipsoid between two circles with authalic latitudes
* &xi;<sub>1</sub> and &xi;<sub>2</sub> is equal to (sin
* &xi;<sub>2</sub> - sin &xi;<sub>1</sub>) \e A / 2, where \e A
* = Area(). For a sphere &xi; = &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &xi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real AuthalicLatitude(real phi) const;
/**
* @param[in] xi the authalic latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* &xi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &phi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real InverseAuthalicLatitude(real xi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &chi; the conformal latitude (degrees).
*
* The conformal latitude, &chi;, gives the mapping of the ellipsoid to a
* sphere which which is conformal (angles are preserved) and in which the
* equator of the ellipsoid maps to the equator of the sphere. For a
* sphere &chi; = &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &chi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real ConformalLatitude(real phi) const;
/**
* @param[in] chi the conformal latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* &chi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold. The returned value
* &phi; lies in [&minus;90&deg;, 90&deg;].
**********************************************************************/
Math::real InverseConformalLatitude(real chi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &psi; the isometric latitude (degrees).
*
* The isometric latitude gives the mapping of the ellipsoid to a plane
* which which is conformal (angles are preserved) and in which the equator
* of the ellipsoid maps to a straight line of constant scale; this mapping
* defines the Mercator projection. For a sphere &psi; =
* sinh<sup>&minus;1</sup> tan &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the result is
* undefined if this condition does not hold. The value returned for &phi;
* = &plusmn;90&deg; is some (positive or negative) large but finite value,
* such that InverseIsometricLatitude returns the original value of &phi;.
**********************************************************************/
Math::real IsometricLatitude(real phi) const;
/**
* @param[in] psi the isometric latitude (degrees).
* @return &phi; the geographic latitude (degrees).
*
* The returned value &phi; lies in [&minus;90&deg;, 90&deg;]. For a
* sphere &phi; = tan<sup>&minus;1</sup> sinh &psi;.
**********************************************************************/
Math::real InverseIsometricLatitude(real psi) const;
///@}
/** \name Other quantities.
**********************************************************************/
///@{
/**
* @param[in] phi the geographic latitude (degrees).
* @return \e R = \e a cos &beta; the radius of a circle of latitude
* &phi; (meters). \e R (&pi;/180&deg;) gives meters per degree
* longitude measured along a circle of latitude.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold.
**********************************************************************/
Math::real CircleRadius(real phi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return \e Z = \e b sin &beta; the distance of a circle of latitude
* &phi; from the equator measured parallel to the ellipsoid axis
* (meters).
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold.
**********************************************************************/
Math::real CircleHeight(real phi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return \e s the distance along a meridian
* between the equator and a point of latitude &phi; (meters). \e s is
* given by \e s = &mu; \e L / 90&deg;, where \e L =
* QuarterMeridian()).
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold.
**********************************************************************/
Math::real MeridianDistance(real phi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &rho; the meridional radius of curvature of the ellipsoid at
* latitude &phi; (meters); this is the curvature of the meridian. \e
* rho is given by &rho; = (180&deg;/&pi;) d\e s / d&phi;,
* where \e s = MeridianDistance(); thus &rho; (&pi;/180&deg;)
* gives meters per degree latitude measured along a meridian.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold.
**********************************************************************/
Math::real MeridionalCurvatureRadius(real phi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @return &nu; the transverse radius of curvature of the ellipsoid at
* latitude &phi; (meters); this is the curvature of a curve on the
* ellipsoid which also lies in a plane perpendicular to the ellipsoid
* and to the meridian. &nu; is related to \e R = CircleRadius() by \e
* R = &nu; cos &phi;.
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the
* result is undefined if this condition does not hold.
**********************************************************************/
Math::real TransverseCurvatureRadius(real phi) const;
/**
* @param[in] phi the geographic latitude (degrees).
* @param[in] azi the angle between the meridian and the normal section
* (degrees).
* @return the radius of curvature of the ellipsoid in the normal
* section at latitude &phi; inclined at an angle \e azi to the
* meridian (meters).
*
* &phi; must lie in the range [&minus;90&deg;, 90&deg;]; the result is
* undefined this condition does not hold.
**********************************************************************/
Math::real NormalCurvatureRadius(real phi, real azi) const;
///@}
/** \name Eccentricity conversions.
**********************************************************************/
///@{
/**
* @param[in] fp = \e f ' = (\e a &minus; \e b) / \e b, the second
* flattening.
* @return \e f = (\e a &minus; \e b) / \e a, the flattening.
*
* \e f ' should lie in (&minus;1, &infin;).
* The returned value \e f lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real SecondFlatteningToFlattening(real fp)
{ return fp / (1 + fp); }
/**
* @param[in] f = (\e a &minus; \e b) / \e a, the flattening.
* @return \e f ' = (\e a &minus; \e b) / \e b, the second flattening.
*
* \e f should lie in (&minus;&infin;, 1).
* The returned value \e f ' lies in (&minus;1, &infin;).
**********************************************************************/
static Math::real FlatteningToSecondFlattening(real f)
{ return f / (1 - f); }
/**
* @param[in] n = (\e a &minus; \e b) / (\e a + \e b), the third
* flattening.
* @return \e f = (\e a &minus; \e b) / \e a, the flattening.
*
* \e n should lie in (&minus;1, 1).
* The returned value \e f lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real ThirdFlatteningToFlattening(real n)
{ return 2 * n / (1 + n); }
/**
* @param[in] f = (\e a &minus; \e b) / \e a, the flattening.
* @return \e n = (\e a &minus; \e b) / (\e a + \e b), the third
* flattening.
*
* \e f should lie in (&minus;&infin;, 1).
* The returned value \e n lies in (&minus;1, 1).
**********************************************************************/
static Math::real FlatteningToThirdFlattening(real f)
{ return f / (2 - f); }
/**
* @param[in] e2 = <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity
* squared.
* @return \e f = (\e a &minus; \e b) / \e a, the flattening.
*
* <i>e</i><sup>2</sup> should lie in (&minus;&infin;, 1).
* The returned value \e f lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real EccentricitySqToFlattening(real e2)
{ using std::sqrt; return e2 / (sqrt(1 - e2) + 1); }
/**
* @param[in] f = (\e a &minus; \e b) / \e a, the flattening.
* @return <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity
* squared.
*
* \e f should lie in (&minus;&infin;, 1).
* The returned value <i>e</i><sup>2</sup> lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real FlatteningToEccentricitySq(real f)
{ return f * (2 - f); }
/**
* @param[in] ep2 = <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity
* squared.
* @return \e f = (\e a &minus; \e b) / \e a, the flattening.
*
* <i>e'</i> <sup>2</sup> should lie in (&minus;1, &infin;).
* The returned value \e f lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real SecondEccentricitySqToFlattening(real ep2)
{ using std::sqrt; return ep2 / (sqrt(1 + ep2) + 1 + ep2); }
/**
* @param[in] f = (\e a &minus; \e b) / \e a, the flattening.
* @return <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity
* squared.
*
* \e f should lie in (&minus;&infin;, 1).
* The returned value <i>e'</i> <sup>2</sup> lies in (&minus;1, &infin;).
**********************************************************************/
static Math::real FlatteningToSecondEccentricitySq(real f)
{ return f * (2 - f) / Math::sq(1 - f); }
/**
* @param[in] epp2 = <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup>
* &minus; <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> +
* <i>b</i><sup>2</sup>), the third eccentricity squared.
* @return \e f = (\e a &minus; \e b) / \e a, the flattening.
*
* <i>e''</i> <sup>2</sup> should lie in (&minus;1, 1).
* The returned value \e f lies in (&minus;&infin;, 1).
**********************************************************************/
static Math::real ThirdEccentricitySqToFlattening(real epp2) {
using std::sqrt;
return 2 * epp2 / (sqrt((1 - epp2) * (1 + epp2)) + 1 + epp2);
}
/**
* @param[in] f = (\e a &minus; \e b) / \e a, the flattening.
* @return <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup> &minus;
* <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> + <i>b</i><sup>2</sup>),
* the third eccentricity squared.
*
* \e f should lie in (&minus;&infin;, 1).
* The returned value <i>e''</i> <sup>2</sup> lies in (&minus;1, 1).
**********************************************************************/
static Math::real FlatteningToThirdEccentricitySq(real f)
{ return f * (2 - f) / (1 + Math::sq(1 - f)); }
///@}
/**
* A global instantiation of Ellipsoid with the parameters for the WGS84
* ellipsoid.
**********************************************************************/
static const Ellipsoid& WGS84();
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_ELLIPSOID_HPP

View File

@@ -0,0 +1,702 @@
/**
* \file EllipticFunction.hpp
* \brief Header for GeographicLib::EllipticFunction class
*
* Copyright (c) Charles Karney (2008-2021) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP)
#define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP 1
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief Elliptic integrals and functions
*
* This provides the elliptic functions and integrals needed for Ellipsoid,
* GeodesicExact, and TransverseMercatorExact. Two categories of function
* are provided:
* - \e static functions to compute symmetric elliptic integrals
* (https://dlmf.nist.gov/19.16.i)
* - \e member functions to compute Legrendre's elliptic
* integrals (https://dlmf.nist.gov/19.2.ii) and the
* Jacobi elliptic functions (https://dlmf.nist.gov/22.2).
* .
* In the latter case, an object is constructed giving the modulus \e k (and
* optionally the parameter &alpha;<sup>2</sup>). The modulus is always
* passed as its square <i>k</i><sup>2</sup> which allows \e k to be pure
* imaginary (<i>k</i><sup>2</sup> &lt; 0). (Confusingly, Abramowitz and
* Stegun call \e m = <i>k</i><sup>2</sup> the "parameter" and \e n =
* &alpha;<sup>2</sup> the "characteristic".)
*
* In geodesic applications, it is convenient to separate the incomplete
* integrals into secular and periodic components, e.g.,
* \f[
* E(\phi, k) = (2 E(k) / \pi) [ \phi + \delta E(\phi, k) ]
* \f]
* where &delta;\e E(&phi;, \e k) is an odd periodic function with period
* &pi;.
*
* The computation of the elliptic integrals uses the algorithms given in
* - B. C. Carlson,
* <a href="https://doi.org/10.1007/BF02198293"> Computation of real or
* complex elliptic integrals</a>, Numerical Algorithms 10, 13--26 (1995)
* .
* with the additional optimizations given in https://dlmf.nist.gov/19.36.i.
* The computation of the Jacobi elliptic functions uses the algorithm given
* in
* - R. Bulirsch,
* <a href="https://doi.org/10.1007/BF01397975"> Numerical Calculation of
* Elliptic Integrals and Elliptic Functions</a>, Numericshe Mathematik 7,
* 78--90 (1965).
* .
* The notation follows https://dlmf.nist.gov/19 and https://dlmf.nist.gov/22
*
* Example of use:
* \include example-EllipticFunction.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT EllipticFunction {
private:
typedef Math::real real;
enum { num_ = 13 }; // Max depth required for sncndn; probably 5 is enough.
real _k2, _kp2, _alpha2, _alphap2, _eps;
real _Kc, _Ec, _Dc, _Pic, _Gc, _Hc;
public:
/** \name Constructor
**********************************************************************/
///@{
/**
* Constructor specifying the modulus and parameter.
*
* @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
* <i>k</i><sup>2</sup> must lie in (&minus;&infin;, 1].
* @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
* &alpha;<sup>2</sup> must lie in (&minus;&infin;, 1].
* @exception GeographicErr if \e k2 or \e alpha2 is out of its legal
* range.
*
* If only elliptic integrals of the first and second kinds are needed,
* then set &alpha;<sup>2</sup> = 0 (the default value); in this case, we
* have &Pi;(&phi;, 0, \e k) = \e F(&phi;, \e k), \e G(&phi;, 0, \e k) = \e
* E(&phi;, \e k), and \e H(&phi;, 0, \e k) = \e F(&phi;, \e k) - \e
* D(&phi;, \e k).
**********************************************************************/
EllipticFunction(real k2 = 0, real alpha2 = 0)
{ Reset(k2, alpha2); }
/**
* Constructor specifying the modulus and parameter and their complements.
*
* @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
* <i>k</i><sup>2</sup> must lie in (&minus;&infin;, 1].
* @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
* &alpha;<sup>2</sup> must lie in (&minus;&infin;, 1].
* @param[in] kp2 the complementary modulus squared <i>k'</i><sup>2</sup> =
* 1 &minus; <i>k</i><sup>2</sup>. This must lie in [0, &infin;).
* @param[in] alphap2 the complementary parameter &alpha;'<sup>2</sup> = 1
* &minus; &alpha;<sup>2</sup>. This must lie in [0, &infin;).
* @exception GeographicErr if \e k2, \e alpha2, \e kp2, or \e alphap2 is
* out of its legal range.
*
* The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2
* = 1. (No checking is done that these conditions are met.) This
* constructor is provided to enable accuracy to be maintained, e.g., when
* \e k is very close to unity.
**********************************************************************/
EllipticFunction(real k2, real alpha2, real kp2, real alphap2)
{ Reset(k2, alpha2, kp2, alphap2); }
/**
* Reset the modulus and parameter.
*
* @param[in] k2 the new value of square of the modulus
* <i>k</i><sup>2</sup> which must lie in (&minus;&infin;, ].
* done.)
* @param[in] alpha2 the new value of parameter &alpha;<sup>2</sup>.
* &alpha;<sup>2</sup> must lie in (&minus;&infin;, 1].
* @exception GeographicErr if \e k2 or \e alpha2 is out of its legal
* range.
**********************************************************************/
void Reset(real k2 = 0, real alpha2 = 0)
{ Reset(k2, alpha2, 1 - k2, 1 - alpha2); }
/**
* Reset the modulus and parameter supplying also their complements.
*
* @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
* <i>k</i><sup>2</sup> must lie in (&minus;&infin;, 1].
* @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
* &alpha;<sup>2</sup> must lie in (&minus;&infin;, 1].
* @param[in] kp2 the complementary modulus squared <i>k'</i><sup>2</sup> =
* 1 &minus; <i>k</i><sup>2</sup>. This must lie in [0, &infin;).
* @param[in] alphap2 the complementary parameter &alpha;'<sup>2</sup> = 1
* &minus; &alpha;<sup>2</sup>. This must lie in [0, &infin;).
* @exception GeographicErr if \e k2, \e alpha2, \e kp2, or \e alphap2 is
* out of its legal range.
*
* The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2
* = 1. (No checking is done that these conditions are met.) This
* constructor is provided to enable accuracy to be maintained, e.g., when
* is very small.
**********************************************************************/
void Reset(real k2, real alpha2, real kp2, real alphap2);
///@}
/** \name Inspector functions.
**********************************************************************/
///@{
/**
* @return the square of the modulus <i>k</i><sup>2</sup>.
**********************************************************************/
Math::real k2() const { return _k2; }
/**
* @return the square of the complementary modulus <i>k'</i><sup>2</sup> =
* 1 &minus; <i>k</i><sup>2</sup>.
**********************************************************************/
Math::real kp2() const { return _kp2; }
/**
* @return the parameter &alpha;<sup>2</sup>.
**********************************************************************/
Math::real alpha2() const { return _alpha2; }
/**
* @return the complementary parameter &alpha;'<sup>2</sup> = 1 &minus;
* &alpha;<sup>2</sup>.
**********************************************************************/
Math::real alphap2() const { return _alphap2; }
///@}
/** \name Complete elliptic integrals.
**********************************************************************/
///@{
/**
* The complete integral of the first kind.
*
* @return \e K(\e k).
*
* \e K(\e k) is defined in https://dlmf.nist.gov/19.2.E4
* \f[
* K(k) = \int_0^{\pi/2} \frac1{\sqrt{1-k^2\sin^2\phi}}\,d\phi.
* \f]
**********************************************************************/
Math::real K() const { return _Kc; }
/**
* The complete integral of the second kind.
*
* @return \e E(\e k).
*
* \e E(\e k) is defined in https://dlmf.nist.gov/19.2.E5
* \f[
* E(k) = \int_0^{\pi/2} \sqrt{1-k^2\sin^2\phi}\,d\phi.
* \f]
**********************************************************************/
Math::real E() const { return _Ec; }
/**
* Jahnke's complete integral.
*
* @return \e D(\e k).
*
* \e D(\e k) is defined in https://dlmf.nist.gov/19.2.E6
* \f[
* D(k) =
* \int_0^{\pi/2} \frac{\sin^2\phi}{\sqrt{1-k^2\sin^2\phi}}\,d\phi.
* \f]
**********************************************************************/
Math::real D() const { return _Dc; }
/**
* The difference between the complete integrals of the first and second
* kinds.
*
* @return \e K(\e k) &minus; \e E(\e k).
**********************************************************************/
Math::real KE() const { return _k2 * _Dc; }
/**
* The complete integral of the third kind.
*
* @return &Pi;(&alpha;<sup>2</sup>, \e k).
*
* &Pi;(&alpha;<sup>2</sup>, \e k) is defined in
* https://dlmf.nist.gov/19.2.E7
* \f[
* \Pi(\alpha^2, k) = \int_0^{\pi/2}
* \frac1{\sqrt{1-k^2\sin^2\phi}(1 - \alpha^2\sin^2\phi)}\,d\phi.
* \f]
**********************************************************************/
Math::real Pi() const { return _Pic; }
/**
* Legendre's complete geodesic longitude integral.
*
* @return \e G(&alpha;<sup>2</sup>, \e k).
*
* \e G(&alpha;<sup>2</sup>, \e k) is given by
* \f[
* G(\alpha^2, k) = \int_0^{\pi/2}
* \frac{\sqrt{1-k^2\sin^2\phi}}{1 - \alpha^2\sin^2\phi}\,d\phi.
* \f]
**********************************************************************/
Math::real G() const { return _Gc; }
/**
* Cayley's complete geodesic longitude difference integral.
*
* @return \e H(&alpha;<sup>2</sup>, \e k).
*
* \e H(&alpha;<sup>2</sup>, \e k) is given by
* \f[
* H(\alpha^2, k) = \int_0^{\pi/2}
* \frac{\cos^2\phi}{(1-\alpha^2\sin^2\phi)\sqrt{1-k^2\sin^2\phi}}
* \,d\phi.
* \f]
**********************************************************************/
Math::real H() const { return _Hc; }
///@}
/** \name Incomplete elliptic integrals.
**********************************************************************/
///@{
/**
* The incomplete integral of the first kind.
*
* @param[in] phi
* @return \e F(&phi;, \e k).
*
* \e F(&phi;, \e k) is defined in https://dlmf.nist.gov/19.2.E4
* \f[
* F(\phi, k) = \int_0^\phi \frac1{\sqrt{1-k^2\sin^2\theta}}\,d\theta.
* \f]
**********************************************************************/
Math::real F(real phi) const;
/**
* The incomplete integral of the second kind.
*
* @param[in] phi
* @return \e E(&phi;, \e k).
*
* \e E(&phi;, \e k) is defined in https://dlmf.nist.gov/19.2.E5
* \f[
* E(\phi, k) = \int_0^\phi \sqrt{1-k^2\sin^2\theta}\,d\theta.
* \f]
**********************************************************************/
Math::real E(real phi) const;
/**
* The incomplete integral of the second kind with the argument given in
* degrees.
*
* @param[in] ang in <i>degrees</i>.
* @return \e E(&pi; <i>ang</i>/180, \e k).
**********************************************************************/
Math::real Ed(real ang) const;
/**
* The inverse of the incomplete integral of the second kind.
*
* @param[in] x
* @return &phi; = <i>E</i><sup>&minus;1</sup>(\e x, \e k); i.e., the
* solution of such that \e E(&phi;, \e k) = \e x.
**********************************************************************/
Math::real Einv(real x) const;
/**
* The incomplete integral of the third kind.
*
* @param[in] phi
* @return &Pi;(&phi;, &alpha;<sup>2</sup>, \e k).
*
* &Pi;(&phi;, &alpha;<sup>2</sup>, \e k) is defined in
* https://dlmf.nist.gov/19.2.E7
* \f[
* \Pi(\phi, \alpha^2, k) = \int_0^\phi
* \frac1{\sqrt{1-k^2\sin^2\theta}(1 - \alpha^2\sin^2\theta)}\,d\theta.
* \f]
**********************************************************************/
Math::real Pi(real phi) const;
/**
* Jahnke's incomplete elliptic integral.
*
* @param[in] phi
* @return \e D(&phi;, \e k).
*
* \e D(&phi;, \e k) is defined in https://dlmf.nist.gov/19.2.E4
* \f[
* D(\phi, k) = \int_0^\phi
* \frac{\sin^2\theta}{\sqrt{1-k^2\sin^2\theta}}\,d\theta.
* \f]
**********************************************************************/
Math::real D(real phi) const;
/**
* Legendre's geodesic longitude integral.
*
* @param[in] phi
* @return \e G(&phi;, &alpha;<sup>2</sup>, \e k).
*
* \e G(&phi;, &alpha;<sup>2</sup>, \e k) is defined by
* \f[
* \begin{align}
* G(\phi, \alpha^2, k) &=
* \frac{k^2}{\alpha^2} F(\phi, k) +
* \biggl(1 - \frac{k^2}{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\
* &= \int_0^\phi
* \frac{\sqrt{1-k^2\sin^2\theta}}{1 - \alpha^2\sin^2\theta}\,d\theta.
* \end{align}
* \f]
*
* Legendre expresses the longitude of a point on the geodesic in terms of
* this combination of elliptic integrals in Exercices de Calcul
* Int&eacute;gral, Vol. 1 (1811), p. 181,
* https://books.google.com/books?id=riIOAAAAQAAJ&pg=PA181.
*
* See \ref geodellip for the expression for the longitude in terms of this
* function.
**********************************************************************/
Math::real G(real phi) const;
/**
* Cayley's geodesic longitude difference integral.
*
* @param[in] phi
* @return \e H(&phi;, &alpha;<sup>2</sup>, \e k).
*
* \e H(&phi;, &alpha;<sup>2</sup>, \e k) is defined by
* \f[
* \begin{align}
* H(\phi, \alpha^2, k) &=
* \frac1{\alpha^2} F(\phi, k) +
* \biggl(1 - \frac1{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\
* &= \int_0^\phi
* \frac{\cos^2\theta}
* {(1-\alpha^2\sin^2\theta)\sqrt{1-k^2\sin^2\theta}}
* \,d\theta.
* \end{align}
* \f]
*
* Cayley expresses the longitude difference of a point on the geodesic in
* terms of this combination of elliptic integrals in Phil. Mag. <b>40</b>
* (1870), p. 333, https://books.google.com/books?id=Zk0wAAAAIAAJ&pg=PA333.
*
* See \ref geodellip for the expression for the longitude in terms of this
* function.
**********************************************************************/
Math::real H(real phi) const;
///@}
/** \name Incomplete integrals in terms of Jacobi elliptic functions.
**********************************************************************/
///@{
/**
* The incomplete integral of the first kind in terms of Jacobi elliptic
* functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return \e F(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
**********************************************************************/
Math::real F(real sn, real cn, real dn) const;
/**
* The incomplete integral of the second kind in terms of Jacobi elliptic
* functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return \e E(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
**********************************************************************/
Math::real E(real sn, real cn, real dn) const;
/**
* The incomplete integral of the third kind in terms of Jacobi elliptic
* functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return &Pi;(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
* (&minus;&pi;, &pi;].
**********************************************************************/
Math::real Pi(real sn, real cn, real dn) const;
/**
* Jahnke's incomplete elliptic integral in terms of Jacobi elliptic
* functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return \e D(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
**********************************************************************/
Math::real D(real sn, real cn, real dn) const;
/**
* Legendre's geodesic longitude integral in terms of Jacobi elliptic
* functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return \e G(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
* (&minus;&pi;, &pi;].
**********************************************************************/
Math::real G(real sn, real cn, real dn) const;
/**
* Cayley's geodesic longitude difference integral in terms of Jacobi
* elliptic functions.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return \e H(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
* (&minus;&pi;, &pi;].
**********************************************************************/
Math::real H(real sn, real cn, real dn) const;
///@}
/** \name Periodic versions of incomplete elliptic integrals.
**********************************************************************/
///@{
/**
* The periodic incomplete integral of the first kind.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; \e F(&phi;, \e k) / (2 \e K(\e k)) -
* &phi;.
**********************************************************************/
Math::real deltaF(real sn, real cn, real dn) const;
/**
* The periodic incomplete integral of the second kind.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; \e E(&phi;, \e k) / (2 \e E(\e k)) -
* &phi;.
**********************************************************************/
Math::real deltaE(real sn, real cn, real dn) const;
/**
* The periodic inverse of the incomplete integral of the second kind.
*
* @param[in] stau = sin&tau;.
* @param[in] ctau = sin&tau;.
* @return the periodic function <i>E</i><sup>&minus;1</sup>(&tau; (2 \e
* E(\e k)/&pi;), \e k) - &tau;.
**********************************************************************/
Math::real deltaEinv(real stau, real ctau) const;
/**
* The periodic incomplete integral of the third kind.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; &Pi;(&phi;, &alpha;<sup>2</sup>,
* \e k) / (2 &Pi;(&alpha;<sup>2</sup>, \e k)) - &phi;.
**********************************************************************/
Math::real deltaPi(real sn, real cn, real dn) const;
/**
* The periodic Jahnke's incomplete elliptic integral.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; \e D(&phi;, \e k) / (2 \e D(\e k)) -
* &phi;.
**********************************************************************/
Math::real deltaD(real sn, real cn, real dn) const;
/**
* Legendre's periodic geodesic longitude integral.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; \e G(&phi;, \e k) / (2 \e G(\e k)) -
* &phi;.
**********************************************************************/
Math::real deltaG(real sn, real cn, real dn) const;
/**
* Cayley's periodic geodesic longitude difference integral.
*
* @param[in] sn = sin&phi;.
* @param[in] cn = cos&phi;.
* @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
* @return the periodic function &pi; \e H(&phi;, \e k) / (2 \e H(\e k)) -
* &phi;.
**********************************************************************/
Math::real deltaH(real sn, real cn, real dn) const;
///@}
/** \name Elliptic functions.
**********************************************************************/
///@{
/**
* The Jacobi elliptic functions.
*
* @param[in] x the argument.
* @param[out] sn sn(\e x, \e k).
* @param[out] cn cn(\e x, \e k).
* @param[out] dn dn(\e x, \e k).
**********************************************************************/
void sncndn(real x, real& sn, real& cn, real& dn) const;
/**
* The &Delta; amplitude function.
*
* @param[in] sn sin&phi;.
* @param[in] cn cos&phi;.
* @return &Delta; = sqrt(1 &minus; <i>k</i><sup>2</sup>
* sin<sup>2</sup>&phi;).
**********************************************************************/
Math::real Delta(real sn, real cn) const {
using std::sqrt;
return sqrt(_k2 < 0 ? 1 - _k2 * sn*sn : _kp2 + _k2 * cn*cn);
}
///@}
/** \name Symmetric elliptic integrals.
**********************************************************************/
///@{
/**
* Symmetric integral of the first kind <i>R</i><sub><i>F</i></sub>.
*
* @param[in] x
* @param[in] y
* @param[in] z
* @return <i>R</i><sub><i>F</i></sub>(\e x, \e y, \e z).
*
* <i>R</i><sub><i>F</i></sub> is defined in https://dlmf.nist.gov/19.16.E1
* \f[ R_F(x, y, z) = \frac12
* \int_0^\infty\frac1{\sqrt{(t + x) (t + y) (t + z)}}\, dt \f]
* If one of the arguments is zero, it is more efficient to call the
* two-argument version of this function with the non-zero arguments.
**********************************************************************/
static real RF(real x, real y, real z);
/**
* Complete symmetric integral of the first kind,
* <i>R</i><sub><i>F</i></sub> with one argument zero.
*
* @param[in] x
* @param[in] y
* @return <i>R</i><sub><i>F</i></sub>(\e x, \e y, 0).
**********************************************************************/
static real RF(real x, real y);
/**
* Degenerate symmetric integral of the first kind
* <i>R</i><sub><i>C</i></sub>.
*
* @param[in] x
* @param[in] y
* @return <i>R</i><sub><i>C</i></sub>(\e x, \e y) =
* <i>R</i><sub><i>F</i></sub>(\e x, \e y, \e y).
*
* <i>R</i><sub><i>C</i></sub> is defined in https://dlmf.nist.gov/19.2.E17
* \f[ R_C(x, y) = \frac12
* \int_0^\infty\frac1{\sqrt{t + x}(t + y)}\,dt \f]
**********************************************************************/
static real RC(real x, real y);
/**
* Symmetric integral of the second kind <i>R</i><sub><i>G</i></sub>.
*
* @param[in] x
* @param[in] y
* @param[in] z
* @return <i>R</i><sub><i>G</i></sub>(\e x, \e y, \e z).
*
* <i>R</i><sub><i>G</i></sub> is defined in Carlson, eq 1.5
* \f[ R_G(x, y, z) = \frac14
* \int_0^\infty[(t + x) (t + y) (t + z)]^{-1/2}
* \biggl(
* \frac x{t + x} + \frac y{t + y} + \frac z{t + z}
* \biggr)t\,dt \f]
* See also https://dlmf.nist.gov/19.16.E3.
* If one of the arguments is zero, it is more efficient to call the
* two-argument version of this function with the non-zero arguments.
**********************************************************************/
static real RG(real x, real y, real z);
/**
* Complete symmetric integral of the second kind,
* <i>R</i><sub><i>G</i></sub> with one argument zero.
*
* @param[in] x
* @param[in] y
* @return <i>R</i><sub><i>G</i></sub>(\e x, \e y, 0).
**********************************************************************/
static real RG(real x, real y);
/**
* Symmetric integral of the third kind <i>R</i><sub><i>J</i></sub>.
*
* @param[in] x
* @param[in] y
* @param[in] z
* @param[in] p
* @return <i>R</i><sub><i>J</i></sub>(\e x, \e y, \e z, \e p).
*
* <i>R</i><sub><i>J</i></sub> is defined in https://dlmf.nist.gov/19.16.E2
* \f[ R_J(x, y, z, p) = \frac32
* \int_0^\infty
* [(t + x) (t + y) (t + z)]^{-1/2} (t + p)^{-1}\, dt \f]
**********************************************************************/
static real RJ(real x, real y, real z, real p);
/**
* Degenerate symmetric integral of the third kind
* <i>R</i><sub><i>D</i></sub>.
*
* @param[in] x
* @param[in] y
* @param[in] z
* @return <i>R</i><sub><i>D</i></sub>(\e x, \e y, \e z) =
* <i>R</i><sub><i>J</i></sub>(\e x, \e y, \e z, \e z).
*
* <i>R</i><sub><i>D</i></sub> is defined in https://dlmf.nist.gov/19.16.E5
* \f[ R_D(x, y, z) = \frac32
* \int_0^\infty[(t + x) (t + y)]^{-1/2} (t + z)^{-3/2}\, dt \f]
**********************************************************************/
static real RD(real x, real y, real z);
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP

143
external/include/GeographicLib/GARS.hpp vendored Normal file
View File

@@ -0,0 +1,143 @@
/**
* \file GARS.hpp
* \brief Header for GeographicLib::GARS class
*
* Copyright (c) Charles Karney (2015-2021) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GARS_HPP)
#define GEOGRAPHICLIB_GARS_HPP 1
#include <GeographicLib/Constants.hpp>
#if defined(_MSC_VER)
// Squelch warnings about dll vs string
# pragma warning (push)
# pragma warning (disable: 4251)
#endif
namespace GeographicLib {
/**
* \brief Conversions for the Global Area Reference System (GARS)
*
* The Global Area Reference System is described in
* - https://en.wikipedia.org/wiki/Global_Area_Reference_System
* - https://earth-info.nga.mil/index.php?dir=coordsys&action=coordsys#tab_gars
* .
* It provides a compact string representation of a geographic area
* (expressed as latitude and longitude). The classes Georef and Geohash
* implement similar compact representations.
*
* Example of use:
* \include example-GARS.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT GARS {
private:
typedef Math::real real;
static const char* const digits_;
static const char* const letters_;
enum {
lonorig_ = -180, // Origin for longitude
latorig_ = -90, // Origin for latitude
baselon_ = 10, // Base for longitude tiles
baselat_ = 24, // Base for latitude tiles
lonlen_ = 3,
latlen_ = 2,
baselen_ = lonlen_ + latlen_,
mult1_ = 2, // base precision = 1/2 degree
mult2_ = 2, // 6th char gives 2x more precision
mult3_ = 3, // 7th char gives 3x more precision
m_ = mult1_ * mult2_ * mult3_,
maxprec_ = 2,
maxlen_ = baselen_ + maxprec_,
};
GARS(); // Disable constructor
public:
/**
* Convert from geographic coordinates to GARS.
*
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[in] prec the precision of the resulting GARS.
* @param[out] gars the GARS string.
* @exception GeographicErr if \e lat is not in [&minus;90&deg;,
* 90&deg;].
* @exception std::bad_alloc if memory for \e gars can't be allocated.
*
* \e prec specifies the precision of \e gars as follows:
* - \e prec = 0 (min), 30' precision, e.g., 006AG;
* - \e prec = 1, 15' precision, e.g., 006AG3;
* - \e prec = 2 (max), 5' precision, e.g., 006AG39.
*
* If \e lat or \e lon is NaN, then \e gars is set to "INVALID".
**********************************************************************/
static void Forward(real lat, real lon, int prec, std::string& gars);
/**
* Convert from GARS to geographic coordinates.
*
* @param[in] gars the GARS.
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] prec the precision of \e gars.
* @param[in] centerp if true (the default) return the center of the
* \e gars, otherwise return the south-west corner.
* @exception GeographicErr if \e gars is illegal.
*
* The case of the letters in \e gars is ignored. \e prec is in the range
* [0, 2] and gives the precision of \e gars as follows:
* - \e prec = 0 (min), 30' precision, e.g., 006AG;
* - \e prec = 1, 15' precision, e.g., 006AG3;
* - \e prec = 2 (max), 5' precision, e.g., 006AG39.
*
* If the first 3 characters of \e gars are "INV", then \e lat and \e lon
* are set to NaN and \e prec is unchanged.
**********************************************************************/
static void Reverse(const std::string& gars, real& lat, real& lon,
int& prec, bool centerp = true);
/**
* The angular resolution of a GARS.
*
* @param[in] prec the precision of the GARS.
* @return the latitude-longitude resolution (degrees).
*
* Internally, \e prec is first put in the range [0, 2].
**********************************************************************/
static Math::real Resolution(int prec) {
return 1/real(prec <= 0 ? mult1_ : (prec == 1 ? mult1_ * mult2_ :
mult1_ * mult2_ * mult3_));
}
/**
* The GARS precision required to meet a given geographic resolution.
*
* @param[in] res the minimum of resolution in latitude and longitude
* (degrees).
* @return GARS precision.
*
* The returned length is in the range [0, 2].
**********************************************************************/
static int Precision(real res) {
using std::abs; res = abs(res);
for (int prec = 0; prec < maxprec_; ++prec)
if (Resolution(prec) <= res)
return prec;
return maxprec_;
}
};
} // namespace GeographicLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // GEOGRAPHICLIB_GARS_HPP

View File

@@ -0,0 +1,553 @@
/**
* \file GeoCoords.hpp
* \brief Header for GeographicLib::GeoCoords class
*
* Copyright (c) Charles Karney (2008-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GEOCOORDS_HPP)
#define GEOGRAPHICLIB_GEOCOORDS_HPP 1
#include <GeographicLib/UTMUPS.hpp>
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief Conversion between geographic coordinates
*
* This class stores a geographic position which may be set via the
* constructors or Reset via
* - latitude and longitude
* - UTM or UPS coordinates
* - a string representation of these or an MGRS coordinate string
*
* The state consists of the latitude and longitude and the supplied UTM or
* UPS coordinates (possibly derived from the MGRS coordinates). If latitude
* and longitude were given then the UTM/UPS coordinates follows the standard
* conventions.
*
* The mutable state consists of the UTM or UPS coordinates for a alternate
* zone. A method SetAltZone is provided to set the alternate UPS/UTM zone.
*
* Methods are provided to return the geographic coordinates, the input UTM
* or UPS coordinates (and associated meridian convergence and scale), or
* alternate UTM or UPS coordinates (and their associated meridian
* convergence and scale).
*
* Once the input string has been parsed, you can print the result out in any
* of the formats, decimal degrees, degrees minutes seconds, MGRS, UTM/UPS.
*
* Example of use:
* \include example-GeoCoords.cpp
*
* <a href="GeoConvert.1.html">GeoConvert</a> is a command-line utility
* providing access to the functionality of GeoCoords.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT GeoCoords {
private:
typedef Math::real real;
real _lat, _long, _easting, _northing, _gamma, _k;
bool _northp;
int _zone; // See UTMUPS::zonespec
mutable real _alt_easting, _alt_northing, _alt_gamma, _alt_k;
mutable int _alt_zone;
void CopyToAlt() const {
_alt_easting = _easting;
_alt_northing = _northing;
_alt_gamma = _gamma;
_alt_k = _k;
_alt_zone = _zone;
}
static void UTMUPSString(int zone, bool northp,
real easting, real northing,
int prec, bool abbrev, std::string& utm);
void FixHemisphere();
public:
/** \name Initializing the GeoCoords object
**********************************************************************/
///@{
/**
* The default constructor sets the coordinate as undefined.
**********************************************************************/
GeoCoords()
: _lat(Math::NaN())
, _long(Math::NaN())
, _easting(Math::NaN())
, _northing(Math::NaN())
, _gamma(Math::NaN())
, _k(Math::NaN())
, _northp(false)
, _zone(UTMUPS::INVALID)
{ CopyToAlt(); }
/**
* Construct from a string.
*
* @param[in] s 1-element, 2-element, or 3-element string representation of
* the position.
* @param[in] centerp governs the interpretation of MGRS coordinates (see
* below).
* @param[in] longfirst governs the interpretation of geographic
* coordinates (see below).
* @exception GeographicErr if the \e s is malformed (see below).
*
* Parse as a string and interpret it as a geographic position. The input
* string is broken into space (or comma) separated pieces and Basic
* decision on which format is based on number of components
* -# MGRS
* -# "Lat Long" or "Long Lat"
* -# "Zone Easting Northing" or "Easting Northing Zone"
*
* The following inputs are approximately the same (Ar Ramadi Bridge, Iraq)
* - Latitude and Longitude
* - 33.44 43.27
* - N33d26.4' E43d16.2'
* - 43d16'12&quot;E 33d26'24&quot;N
* - 43:16:12E 33:26:24
* - MGRS
* - 38SLC30
* - 38SLC391014
* - 38SLC3918701405
* - 37SHT9708
* - UTM
* - 38n 339188 3701405
* - 897039 3708229 37n
*
* <b>Latitude and Longitude parsing</b>: Latitude precedes longitude,
* unless a N, S, E, W hemisphere designator is used on one or both
* coordinates. If \e longfirst = true (default is false), then
* longitude precedes latitude in the absence of a hemisphere designator.
* Thus (with \e longfirst = false)
* - 40 -75
* - N40 W75
* - -75 N40
* - 75W 40N
* - E-75 -40S
* .
* are all the same position. The coordinates may be given in
* decimal degrees, degrees and decimal minutes, degrees, minutes,
* seconds, etc. Use d, ', and &quot; to mark off the degrees,
* minutes and seconds. Various alternative symbols for degrees, minutes,
* and seconds are allowed. Alternatively, use : to separate these
* components. A single addition or subtraction is allowed. (See
* DMS::Decode for details.) Thus
* - 40d30'30&quot;
* - 40d30'30
* - 40&deg;30'30
* - 40d30.5'
* - 40d30.5
* - 40:30:30
* - 40:30.5
* - 40.508333333
* - 40:30+0:0:30
* - 40:31-0:0.5
* .
* all specify the same angle. The leading sign applies to the following
* components so -1d30 is -(1+30/60) = &minus;1.5. However, note
* that -1:30-0:0:15 is parsed as (-1:30) + (-0:0:15) = &minus;(1+30/60)
* &minus; (15/3600). Latitudes must be in the range [&minus;90&deg;,
* 90&deg;]. Internally longitudes are reduced to the range
* [&minus;180&deg;, 180&deg;].
*
* <b>UTM/UPS parsing</b>: For UTM zones (&minus;80&deg; &le; Lat <
* 84&deg;), the zone designator is made up of a zone number (for 1 to 60)
* and a hemisphere letter (n or s), e.g., 38n (38north can also be used).
* The latitude band designer ([C--M] in the southern hemisphere and [N--X]
* in the northern) should NOT be used. (This is part of the MGRS
* coordinate.) The zone designator for the poles (where UPS is employed)
* is a hemisphere letter by itself, i.e., n or s (north or south can also
* be used).
*
* <b>MGRS parsing</b> interprets the grid references as square area at the
* specified precision (1m, 10m, 100m, etc.). If \e centerp = true (the
* default), the center of this square is then taken to be the precise
* position; thus:
* - 38SMB = 38n 450000 3650000
* - 38SMB4484 = 38n 444500 3684500
* - 38SMB44148470 = 38n 444145 3684705
* .
* Otherwise, the "south-west" corner of the square is used, i.e.,
* - 38SMB = 38n 400000 3600000
* - 38SMB4484 = 38n 444000 3684000
* - 38SMB44148470 = 38n 444140 3684700
**********************************************************************/
explicit GeoCoords(const std::string& s,
bool centerp = true, bool longfirst = false)
{ Reset(s, centerp, longfirst); }
/**
* Construct from geographic coordinates.
*
* @param[in] latitude (degrees).
* @param[in] longitude (degrees).
* @param[in] zone if specified, force the UTM/UPS representation to use a
* specified zone using the rules given in UTMUPS::zonespec.
* @exception GeographicErr if \e latitude is not in [&minus;90&deg;,
* 90&deg;].
* @exception GeographicErr if \e zone cannot be used for this location.
**********************************************************************/
GeoCoords(real latitude, real longitude, int zone = UTMUPS::STANDARD) {
Reset(latitude, longitude, zone);
}
/**
* Construct from UTM/UPS coordinates.
*
* @param[in] zone UTM zone (zero means UPS).
* @param[in] northp hemisphere (true means north, false means south).
* @param[in] easting (meters).
* @param[in] northing (meters).
* @exception GeographicErr if \e zone, \e easting, or \e northing is
* outside its allowed range.
**********************************************************************/
GeoCoords(int zone, bool northp, real easting, real northing) {
Reset(zone, northp, easting, northing);
}
/**
* Reset the location from a string. See
* GeoCoords(const std::string& s, bool centerp, bool longfirst).
*
* @param[in] s 1-element, 2-element, or 3-element string representation of
* the position.
* @param[in] centerp governs the interpretation of MGRS coordinates.
* @param[in] longfirst governs the interpretation of geographic
* coordinates.
* @exception GeographicErr if the \e s is malformed.
**********************************************************************/
void Reset(const std::string& s,
bool centerp = true, bool longfirst = false);
/**
* Reset the location in terms of geographic coordinates. See
* GeoCoords(real latitude, real longitude, int zone).
*
* @param[in] latitude (degrees).
* @param[in] longitude (degrees).
* @param[in] zone if specified, force the UTM/UPS representation to use a
* specified zone using the rules given in UTMUPS::zonespec.
* @exception GeographicErr if \e latitude is not in [&minus;90&deg;,
* 90&deg;].
* @exception GeographicErr if \e zone cannot be used for this location.
**********************************************************************/
void Reset(real latitude, real longitude, int zone = UTMUPS::STANDARD) {
UTMUPS::Forward(latitude, longitude,
_zone, _northp, _easting, _northing, _gamma, _k,
zone);
_lat = latitude;
_long = longitude;
if (_long >= 180) _long -= 360;
else if (_long < -180) _long += 360;
CopyToAlt();
}
/**
* Reset the location in terms of UPS/UPS coordinates. See
* GeoCoords(int zone, bool northp, real easting, real northing).
*
* @param[in] zone UTM zone (zero means UPS).
* @param[in] northp hemisphere (true means north, false means south).
* @param[in] easting (meters).
* @param[in] northing (meters).
* @exception GeographicErr if \e zone, \e easting, or \e northing is
* outside its allowed range.
**********************************************************************/
void Reset(int zone, bool northp, real easting, real northing) {
UTMUPS::Reverse(zone, northp, easting, northing,
_lat, _long, _gamma, _k);
_zone = zone;
_northp = northp;
_easting = easting;
_northing = northing;
FixHemisphere();
CopyToAlt();
}
///@}
/** \name Querying the GeoCoords object
**********************************************************************/
///@{
/**
* @return latitude (degrees)
**********************************************************************/
Math::real Latitude() const { return _lat; }
/**
* @return longitude (degrees)
**********************************************************************/
Math::real Longitude() const { return _long; }
/**
* @return easting (meters)
**********************************************************************/
Math::real Easting() const { return _easting; }
/**
* @return northing (meters)
**********************************************************************/
Math::real Northing() const { return _northing; }
/**
* @return meridian convergence (degrees) for the UTM/UPS projection.
**********************************************************************/
Math::real Convergence() const { return _gamma; }
/**
* @return scale for the UTM/UPS projection.
**********************************************************************/
Math::real Scale() const { return _k; }
/**
* @return hemisphere (false means south, true means north).
**********************************************************************/
bool Northp() const { return _northp; }
/**
* @return hemisphere letter n or s.
**********************************************************************/
char Hemisphere() const { return _northp ? 'n' : 's'; }
/**
* @return the zone corresponding to the input (return 0 for UPS).
**********************************************************************/
int Zone() const { return _zone; }
///@}
/** \name Setting and querying the alternate zone
**********************************************************************/
///@{
/**
* Specify alternate zone number.
*
* @param[in] zone zone number for the alternate representation.
* @exception GeographicErr if \e zone cannot be used for this location.
*
* See UTMUPS::zonespec for more information on the interpretation of \e
* zone. Note that \e zone == UTMUPS::STANDARD (the default) use the
* standard UPS or UTM zone, UTMUPS::MATCH does nothing retaining the
* existing alternate representation. Before this is called the alternate
* zone is the input zone.
**********************************************************************/
void SetAltZone(int zone = UTMUPS::STANDARD) const {
if (zone == UTMUPS::MATCH)
return;
zone = UTMUPS::StandardZone(_lat, _long, zone);
if (zone == _zone)
CopyToAlt();
else {
bool northp;
UTMUPS::Forward(_lat, _long,
_alt_zone, northp,
_alt_easting, _alt_northing, _alt_gamma, _alt_k,
zone);
}
}
/**
* @return current alternate zone (return 0 for UPS).
**********************************************************************/
int AltZone() const { return _alt_zone; }
/**
* @return easting (meters) for alternate zone.
**********************************************************************/
Math::real AltEasting() const { return _alt_easting; }
/**
* @return northing (meters) for alternate zone.
**********************************************************************/
Math::real AltNorthing() const { return _alt_northing; }
/**
* @return meridian convergence (degrees) for alternate zone.
**********************************************************************/
Math::real AltConvergence() const { return _alt_gamma; }
/**
* @return scale for alternate zone.
**********************************************************************/
Math::real AltScale() const { return _alt_k; }
///@}
/** \name String representations of the GeoCoords object
**********************************************************************/
///@{
/**
* String representation with latitude and longitude as signed decimal
* degrees.
*
* @param[in] prec precision (relative to about 1m).
* @param[in] longfirst if true give longitude first (default = false)
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return decimal latitude/longitude string representation.
*
* Precision specifies accuracy of representation as follows:
* - prec = &minus;5 (min), 1&deg;
* - prec = 0, 10<sup>&minus;5</sup>&deg; (about 1m)
* - prec = 3, 10<sup>&minus;8</sup>&deg;
* - prec = 9 (max), 10<sup>&minus;14</sup>&deg;
**********************************************************************/
std::string GeoRepresentation(int prec = 0, bool longfirst = false) const;
/**
* String representation with latitude and longitude as degrees, minutes,
* seconds, and hemisphere.
*
* @param[in] prec precision (relative to about 1m)
* @param[in] longfirst if true give longitude first (default = false)
* @param[in] dmssep if non-null, use as the DMS separator character
* (instead of d, ', &quot; delimiters).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return DMS latitude/longitude string representation.
*
* Precision specifies accuracy of representation as follows:
* - prec = &minus;5 (min), 1&deg;
* - prec = &minus;4, 0.1&deg;
* - prec = &minus;3, 1'
* - prec = &minus;2, 0.1'
* - prec = &minus;1, 1&quot;
* - prec = 0, 0.1&quot; (about 3m)
* - prec = 1, 0.01&quot;
* - prec = 10 (max), 10<sup>&minus;11</sup>&quot;
**********************************************************************/
std::string DMSRepresentation(int prec = 0, bool longfirst = false,
char dmssep = char(0))
const;
/**
* MGRS string.
*
* @param[in] prec precision (relative to about 1m).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return MGRS string.
*
* This gives the coordinates of the enclosing grid square with size given
* by the precision. Thus 38n 444180 3684790 converted to a MGRS
* coordinate at precision &minus;2 (100m) is 38SMB441847 and not
* 38SMB442848. \e prec specifies the precision of the MGRS string as
* follows:
* - prec = &minus;6 (min), only the grid zone is returned, e.g., 38S
* - prec = &minus;5, 100km, e.g., 38SMB
* - prec = &minus;4, 10km
* - prec = &minus;3, 1km
* - prec = &minus;2, 100m
* - prec = &minus;1, 10m
* - prec = 0, 1m
* - prec = 1, 0.1m
* - prec = 6 (max), 1&mu;m
**********************************************************************/
std::string MGRSRepresentation(int prec = 0) const;
/**
* UTM/UPS string.
*
* @param[in] prec precision (relative to about 1m)
* @param[in] abbrev if true (the default) use abbreviated (n/s) notation
* for hemisphere; otherwise spell out the hemisphere (north/south)
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return UTM/UPS string representation: zone designator, easting, and
* northing.
*
* Precision specifies accuracy of representation as follows:
* - prec = &minus;5 (min), 100km
* - prec = &minus;3, 1km
* - prec = 0, 1m
* - prec = 3, 1mm
* - prec = 6, 1&mu;m
* - prec = 9 (max), 1nm
**********************************************************************/
std::string UTMUPSRepresentation(int prec = 0, bool abbrev = true) const;
/**
* UTM/UPS string with hemisphere override.
*
* @param[in] northp hemisphere override
* @param[in] prec precision (relative to about 1m)
* @param[in] abbrev if true (the default) use abbreviated (n/s) notation
* for hemisphere; otherwise spell out the hemisphere (north/south)
* @exception GeographicErr if the hemisphere override attempts to change
* UPS N to UPS S or vice versa.
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return UTM/UPS string representation: zone designator, easting, and
* northing.
**********************************************************************/
std::string UTMUPSRepresentation(bool northp, int prec = 0,
bool abbrev = true) const;
/**
* MGRS string for the alternate zone. See GeoCoords::MGRSRepresentation.
*
* @param[in] prec precision (relative to about 1m).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return MGRS string.
**********************************************************************/
std::string AltMGRSRepresentation(int prec = 0) const;
/**
* UTM/UPS string for the alternate zone. See
* GeoCoords::UTMUPSRepresentation.
*
* @param[in] prec precision (relative to about 1m)
* @param[in] abbrev if true (the default) use abbreviated (n/s) notation
* for hemisphere; otherwise spell out the hemisphere (north/south)
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return UTM/UPS string representation: zone designator, easting, and
* northing.
**********************************************************************/
std::string AltUTMUPSRepresentation(int prec = 0, bool abbrev = true)
const;
/**
* UTM/UPS string for the alternate zone, with hemisphere override.
*
* @param[in] northp hemisphere override
* @param[in] prec precision (relative to about 1m)
* @param[in] abbrev if true (the default) use abbreviated (n/s) notation
* for hemisphere; otherwise spell out the hemisphere (north/south)
* @exception GeographicErr if the hemisphere override attempts to change
* UPS n to UPS s or vice verse.
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return UTM/UPS string representation: zone designator, easting, and
* northing.
**********************************************************************/
std::string AltUTMUPSRepresentation(bool northp, int prec = 0,
bool abbrev = true) const;
///@}
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
*
* (The WGS84 value is returned because the UTM and UPS projections are
* based on this ellipsoid.)
**********************************************************************/
Math::real EquatorialRadius() const { return UTMUPS::EquatorialRadius(); }
/**
* @return \e f the flattening of the WGS84 ellipsoid.
*
* (The WGS84 value is returned because the UTM and UPS projections are
* based on this ellipsoid.)
**********************************************************************/
Math::real Flattening() const { return UTMUPS::Flattening(); }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_GEOCOORDS_HPP

View File

@@ -0,0 +1,274 @@
/**
* \file Geocentric.hpp
* \brief Header for GeographicLib::Geocentric class
*
* Copyright (c) Charles Karney (2008-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP)
#define GEOGRAPHICLIB_GEOCENTRIC_HPP 1
#include <vector>
#include <GeographicLib/Constants.hpp>
namespace GeographicLib {
/**
* \brief %Geocentric coordinates
*
* Convert between geodetic coordinates latitude = \e lat, longitude = \e
* lon, height = \e h (measured vertically from the surface of the ellipsoid)
* to geocentric coordinates (\e X, \e Y, \e Z). The origin of geocentric
* coordinates is at the center of the earth. The \e Z axis goes thru the
* north pole, \e lat = 90&deg;. The \e X axis goes thru \e lat = 0,
* \e lon = 0. %Geocentric coordinates are also known as earth centered,
* earth fixed (ECEF) coordinates.
*
* The conversion from geographic to geocentric coordinates is
* straightforward. For the reverse transformation we use
* - H. Vermeille,
* <a href="https://doi.org/10.1007/s00190-002-0273-6"> Direct
* transformation from geocentric coordinates to geodetic coordinates</a>,
* J. Geodesy 76, 451--454 (2002).
* .
* Several changes have been made to ensure that the method returns accurate
* results for all finite inputs (even if \e h is infinite). The changes are
* described in Appendix B of
* - C. F. F. Karney,
* <a href="https://arxiv.org/abs/1102.1215v1">Geodesics
* on an ellipsoid of revolution</a>,
* Feb. 2011;
* preprint
* <a href="https://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.
* .
* Vermeille similarly updated his method in
* - H. Vermeille,
* <a href="https://doi.org/10.1007/s00190-010-0419-x">
* An analytical method to transform geocentric into
* geodetic coordinates</a>, J. Geodesy 85, 105--117 (2011).
* .
* See \ref geocentric for more information.
*
* The errors in these routines are close to round-off. Specifically, for
* points within 5000 km of the surface of the ellipsoid (either inside or
* outside the ellipsoid), the error is bounded by 7 nm (7 nanometers) for
* the WGS84 ellipsoid. See \ref geocentric for further information on the
* errors.
*
* Example of use:
* \include example-Geocentric.cpp
*
* <a href="CartConvert.1.html">CartConvert</a> is a command-line utility
* providing access to the functionality of Geocentric and LocalCartesian.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT Geocentric {
private:
typedef Math::real real;
friend class LocalCartesian;
friend class MagneticCircle; // MagneticCircle uses Rotation
friend class MagneticModel; // MagneticModel uses IntForward
friend class GravityCircle; // GravityCircle uses Rotation
friend class GravityModel; // GravityModel uses IntForward
friend class NormalGravity; // NormalGravity uses IntForward
static const size_t dim_ = 3;
static const size_t dim2_ = dim_ * dim_;
real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad;
static void Rotation(real sphi, real cphi, real slam, real clam,
real M[dim2_]);
static void Rotate(const real M[dim2_], real x, real y, real z,
real& X, real& Y, real& Z) {
// Perform [X,Y,Z]^t = M.[x,y,z]^t
// (typically local cartesian to geocentric)
X = M[0] * x + M[1] * y + M[2] * z;
Y = M[3] * x + M[4] * y + M[5] * z;
Z = M[6] * x + M[7] * y + M[8] * z;
}
static void Unrotate(const real M[dim2_], real X, real Y, real Z,
real& x, real& y, real& z) {
// Perform [x,y,z]^t = M^t.[X,Y,Z]^t
// (typically geocentric to local cartesian)
x = M[0] * X + M[3] * Y + M[6] * Z;
y = M[1] * X + M[4] * Y + M[7] * Z;
z = M[2] * X + M[5] * Y + M[8] * Z;
}
void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z,
real M[dim2_]) const;
void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h,
real M[dim2_]) const;
public:
/**
* Constructor for a ellipsoid with
*
* @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
* positive.
**********************************************************************/
Geocentric(real a, real f);
/**
* A default constructor (for use by NormalGravity).
**********************************************************************/
Geocentric() : _a(-1) {}
/**
* Convert from geodetic to geocentric coordinates.
*
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[in] h height of point above the ellipsoid (meters).
* @param[out] X geocentric coordinate (meters).
* @param[out] Y geocentric coordinate (meters).
* @param[out] Z geocentric coordinate (meters).
*
* \e lat should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
void Forward(real lat, real lon, real h, real& X, real& Y, real& Z)
const {
if (Init())
IntForward(lat, lon, h, X, Y, Z, NULL);
}
/**
* Convert from geodetic to geocentric coordinates and return rotation
* matrix.
*
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[in] h height of point above the ellipsoid (meters).
* @param[out] X geocentric coordinate (meters).
* @param[out] Y geocentric coordinate (meters).
* @param[out] Z geocentric coordinate (meters).
* @param[out] M if the length of the vector is 9, fill with the rotation
* matrix in row-major order.
*
* Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
* express \e v as \e column vectors in one of two ways
* - in east, north, up coordinates (where the components are relative to a
* local coordinate system at (\e lat, \e lon, \e h)); call this
* representation \e v1.
* - in geocentric \e X, \e Y, \e Z coordinates; call this representation
* \e v0.
* .
* Then we have \e v0 = \e M &sdot; \e v1.
**********************************************************************/
void Forward(real lat, real lon, real h, real& X, real& Y, real& Z,
std::vector<real>& M)
const {
if (!Init())
return;
if (M.end() == M.begin() + dim2_) {
real t[dim2_];
IntForward(lat, lon, h, X, Y, Z, t);
std::copy(t, t + dim2_, M.begin());
} else
IntForward(lat, lon, h, X, Y, Z, NULL);
}
/**
* Convert from geocentric to geodetic to coordinates.
*
* @param[in] X geocentric coordinate (meters).
* @param[in] Y geocentric coordinate (meters).
* @param[in] Z geocentric coordinate (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters).
*
* In general, there are multiple solutions and the result which minimizes
* |<i>h</i> |is returned, i.e., (<i>lat</i>, <i>lon</i>) corresponds to
* the closest point on the ellipsoid. If there are still multiple
* solutions with different latitudes (applies only if \e Z = 0), then the
* solution with \e lat > 0 is returned. If there are still multiple
* solutions with different longitudes (applies only if \e X = \e Y = 0)
* then \e lon = 0 is returned. The value of \e h returned satisfies \e h
* &ge; &minus; \e a (1 &minus; <i>e</i><sup>2</sup>) / sqrt(1 &minus;
* <i>e</i><sup>2</sup> sin<sup>2</sup>\e lat). The value of \e lon
* returned is in the range [&minus;180&deg;, 180&deg;].
**********************************************************************/
void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h)
const {
if (Init())
IntReverse(X, Y, Z, lat, lon, h, NULL);
}
/**
* Convert from geocentric to geodetic to coordinates.
*
* @param[in] X geocentric coordinate (meters).
* @param[in] Y geocentric coordinate (meters).
* @param[in] Z geocentric coordinate (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters).
* @param[out] M if the length of the vector is 9, fill with the rotation
* matrix in row-major order.
*
* Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
* express \e v as \e column vectors in one of two ways
* - in east, north, up coordinates (where the components are relative to a
* local coordinate system at (\e lat, \e lon, \e h)); call this
* representation \e v1.
* - in geocentric \e X, \e Y, \e Z coordinates; call this representation
* \e v0.
* .
* Then we have \e v1 = <i>M</i><sup>T</sup> &sdot; \e v0, where
* <i>M</i><sup>T</sup> is the transpose of \e M.
**********************************************************************/
void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h,
std::vector<real>& M)
const {
if (!Init())
return;
if (M.end() == M.begin() + dim2_) {
real t[dim2_];
IntReverse(X, Y, Z, lat, lon, h, t);
std::copy(t, t + dim2_, M.begin());
} else
IntReverse(X, Y, Z, lat, lon, h, NULL);
}
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return true if the object has been initialized.
**********************************************************************/
bool Init() const { return _a > 0; }
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const
{ return Init() ? _a : Math::NaN(); }
/**
* @return \e f the flattening of the ellipsoid. This is the
* value used in the constructor.
**********************************************************************/
Math::real Flattening() const
{ return Init() ? _f : Math::NaN(); }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
/**
* A global instantiation of Geocentric with the parameters for the WGS84
* ellipsoid.
**********************************************************************/
static const Geocentric& WGS84();
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_GEOCENTRIC_HPP

View File

@@ -0,0 +1,977 @@
/**
* \file Geodesic.hpp
* \brief Header for GeographicLib::Geodesic class
*
* Copyright (c) Charles Karney (2009-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESIC_HPP)
#define GEOGRAPHICLIB_GEODESIC_HPP 1
#include <GeographicLib/Constants.hpp>
#if !defined(GEOGRAPHICLIB_GEODESIC_ORDER)
/**
* The order of the expansions used by Geodesic.
* GEOGRAPHICLIB_GEODESIC_ORDER can be set to any integer in [3, 8].
**********************************************************************/
# define GEOGRAPHICLIB_GEODESIC_ORDER \
(GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
(GEOGRAPHICLIB_PRECISION == 1 ? 3 : \
(GEOGRAPHICLIB_PRECISION == 3 ? 7 : 8)))
#endif
namespace GeographicLib {
class GeodesicLine;
/**
* \brief %Geodesic calculations
*
* The shortest path between two points on a ellipsoid at (\e lat1, \e lon1)
* and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and
* the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
* the two end points. (The azimuth is the heading measured clockwise from
* north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you
* beyond point 2 not back to point 1.) In the figure below, latitude if
* labeled &phi;, longitude &lambda; (with &lambda;<sub>12</sub> =
* &lambda;<sub>2</sub> &minus; &lambda;<sub>1</sub>), and azimuth &alpha;.
*
* <img src="https://upload.wikimedia.org/wikipedia/commons/c/cb/Geodesic_problem_on_an_ellipsoid.svg" width=250 alt="spheroidal triangle">
*
* Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
* lon2, and \e azi2. This is the \e direct geodesic problem and its
* solution is given by the function Geodesic::Direct. (If \e s12 is
* sufficiently large that the geodesic wraps more than halfway around the
* earth, there will be another geodesic between the points with a smaller \e
* s12.)
*
* Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e
* azi2, and \e s12. This is the \e inverse geodesic problem, whose solution
* is given by Geodesic::Inverse. Usually, the solution to the inverse
* problem is unique. In cases where there are multiple solutions (all with
* the same \e s12, of course), all the solutions can be easily generated
* once a particular solution is provided.
*
* The standard way of specifying the direct problem is the specify the
* distance \e s12 to the second point. However it is sometimes useful
* instead to specify the arc length \e a12 (in degrees) on the auxiliary
* sphere. This is a mathematical construct used in solving the geodesic
* problems. The solution of the direct problem in this form is provided by
* Geodesic::ArcDirect. An arc length in excess of 180&deg; indicates that
* the geodesic is not a shortest path. In addition, the arc length between
* an equatorial crossing and the next extremum of latitude for a geodesic is
* 90&deg;.
*
* This class can also calculate several other quantities related to
* geodesics. These are:
* - <i>reduced length</i>. If we fix the first point and increase \e azi1
* by \e dazi1 (radians), the second point is displaced \e m12 \e dazi1 in
* the direction \e azi2 + 90&deg;. The quantity \e m12 is called
* the "reduced length" and is symmetric under interchange of the two
* points. On a curved surface the reduced length obeys a symmetry
* relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e
* s12. The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an
* azimuthal equidistant projection.
* - <i>geodesic scale</i>. Consider a reference geodesic and a second
* geodesic parallel to this one at point 1 and separated by a small
* distance \e dt. The separation of the two geodesics at point 2 is \e
* M12 \e dt where \e M12 is called the "geodesic scale". \e M21 is
* defined similarly (with the geodesics being parallel at point 2). On a
* flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gives
* the scale of the Cassini-Soldner projection.
* - <i>area</i>. The area between the geodesic from point 1 to point 2 and
* the equation is represented by \e S12; it is the area, measured
* counter-clockwise, of the geodesic quadrilateral with corners
* (<i>lat1</i>,<i>lon1</i>), (0,<i>lon1</i>), (0,<i>lon2</i>), and
* (<i>lat2</i>,<i>lon2</i>). It can be used to compute the area of any
* geodesic polygon.
*
* Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and
* Geodesic::Inverse allow these quantities to be returned. In addition
* there are general functions Geodesic::GenDirect, and Geodesic::GenInverse
* which allow an arbitrary set of results to be computed. The quantities \e
* m12, \e M12, \e M21 which all specify the behavior of nearby geodesics
* obey addition rules. If points 1, 2, and 3 all lie on a single geodesic,
* then the following rules hold:
* - \e s13 = \e s12 + \e s23
* - \e a13 = \e a12 + \e a23
* - \e S13 = \e S12 + \e S23
* - \e m13 = \e m12 \e M23 + \e m23 \e M21
* - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e m23 / \e m12
* - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e m12 / \e m23
*
* Additional functionality is provided by the GeodesicLine class, which
* allows a sequence of points along a geodesic to be computed.
*
* The shortest distance returned by the solution of the inverse problem is
* (obviously) uniquely defined. However, in a few special cases there are
* multiple azimuths which yield the same shortest distance. Here is a
* catalog of those cases:
* - \e lat1 = &minus;\e lat2 (with neither point at a pole). If \e azi1 =
* \e azi2, the geodesic is unique. Otherwise there are two geodesics and
* the second one is obtained by setting [\e azi1, \e azi2] &rarr; [\e
* azi2, \e azi1], [\e M12, \e M21] &rarr; [\e M21, \e M12], \e S12 &rarr;
* &minus;\e S12. (This occurs when the longitude difference is near
* &plusmn;180&deg; for oblate ellipsoids.)
* - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole). If
* \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique. Otherwise
* there are two geodesics and the second one is obtained by setting [\e
* azi1, \e azi2] &rarr; [&minus;\e azi1, &minus;\e azi2], \e S12 &rarr;
* &minus;\e S12. (This occurs when \e lat2 is near &minus;\e lat1 for
* prolate ellipsoids.)
* - Points 1 and 2 at opposite poles. There are infinitely many geodesics
* which can be generated by setting [\e azi1, \e azi2] &rarr; [\e azi1, \e
* azi2] + [\e d, &minus;\e d], for arbitrary \e d. (For spheres, this
* prescription applies when points 1 and 2 are antipodal.)
* - \e s12 = 0 (coincident points). There are infinitely many geodesics
* which can be generated by setting [\e azi1, \e azi2] &rarr;
* [\e azi1, \e azi2] + [\e d, \e d], for arbitrary \e d.
*
* The calculations are accurate to better than 15 nm (15 nanometers) for the
* WGS84 ellipsoid. See Sec. 9 of
* <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for
* details. The algorithms used by this class are based on series expansions
* using the flattening \e f as a small parameter. These are only accurate
* for |<i>f</i>| &lt; 0.02; however reasonably accurate results will be
* obtained for |<i>f</i>| &lt; 0.2. Here is a table of the approximate
* maximum error (expressed as a distance) for an ellipsoid with the same
* equatorial radius as the WGS84 ellipsoid and different values of the
* flattening.<pre>
* |f| error
* 0.01 25 nm
* 0.02 30 nm
* 0.05 10 um
* 0.1 1.5 mm
* 0.2 300 mm
* </pre>
* For very eccentric ellipsoids, use GeodesicExact instead.
*
* The algorithms are described in
* - C. F. F. Karney,
* <a href="https://doi.org/10.1007/s00190-012-0578-z">
* Algorithms for geodesics</a>,
* J. Geodesy <b>87</b>, 43--55 (2013);
* DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
* 10.1007/s00190-012-0578-z</a>;
* addenda:
* <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
* geod-addenda.html</a>.
* .
* For more information on geodesics see \ref geodesic.
*
* Example of use:
* \include example-Geodesic.cpp
*
* <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
* providing access to the functionality of Geodesic and GeodesicLine.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT Geodesic {
private:
typedef Math::real real;
friend class GeodesicLine;
static const int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nA3x_ = nA3_;
static const int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
static const int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
// Size for temporary array
// nC = max(max(nC1_, nC1p_, nC2_) + 1, max(nC3_, nC4_))
static const int nC_ = GEOGRAPHICLIB_GEODESIC_ORDER + 1;
static const unsigned maxit1_ = 20;
unsigned maxit2_;
real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
enum captype {
CAP_NONE = 0U,
CAP_C1 = 1U<<0,
CAP_C1p = 1U<<1,
CAP_C2 = 1U<<2,
CAP_C3 = 1U<<3,
CAP_C4 = 1U<<4,
CAP_ALL = 0x1FU,
CAP_MASK = CAP_ALL,
OUT_ALL = 0x7F80U,
OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
};
static real SinCosSeries(bool sinp,
real sinx, real cosx, const real c[], int n);
static real Astroid(real x, real y);
real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_];
void Lengths(real eps, real sig12,
real ssig1, real csig1, real dn1,
real ssig2, real csig2, real dn2,
real cbet1, real cbet2, unsigned outmask,
real& s12s, real& m12a, real& m0,
real& M12, real& M21, real Ca[]) const;
real InverseStart(real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real lam12, real slam12, real clam12,
real& salp1, real& calp1,
real& salp2, real& calp2, real& dnm,
real Ca[]) const;
real Lambda12(real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real salp1, real calp1, real slam120, real clam120,
real& salp2, real& calp2, real& sig12,
real& ssig1, real& csig1, real& ssig2, real& csig2,
real& eps, real& domg12,
bool diffp, real& dlam12, real Ca[]) const;
real GenInverse(real lat1, real lon1, real lat2, real lon2,
unsigned outmask, real& s12,
real& salp1, real& calp1, real& salp2, real& calp2,
real& m12, real& M12, real& M21, real& S12) const;
// These are Maxima generated functions to provide series approximations to
// the integrals for the ellipsoidal geodesic.
static real A1m1f(real eps);
static void C1f(real eps, real c[]);
static void C1pf(real eps, real c[]);
static real A2m1f(real eps);
static void C2f(real eps, real c[]);
void A3coeff();
real A3f(real eps) const;
void C3coeff();
void C3f(real eps, real c[]) const;
void C4coeff();
void C4f(real k2, real c[]) const;
public:
/**
* Bit masks for what calculations to do. These masks do double duty.
* They signify to the GeodesicLine::GeodesicLine constructor and to
* Geodesic::Line what capabilities should be included in the GeodesicLine
* object. They also specify which results to return in the general
* routines Geodesic::GenDirect and Geodesic::GenInverse routines.
* GeodesicLine::mask is a duplication of this enum.
**********************************************************************/
enum mask {
/**
* No capabilities, no output.
* @hideinitializer
**********************************************************************/
NONE = 0U,
/**
* Calculate latitude \e lat2. (It's not necessary to include this as a
* capability to GeodesicLine because this is included by default.)
* @hideinitializer
**********************************************************************/
LATITUDE = 1U<<7 | CAP_NONE,
/**
* Calculate longitude \e lon2.
* @hideinitializer
**********************************************************************/
LONGITUDE = 1U<<8 | CAP_C3,
/**
* Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
* include this as a capability to GeodesicLine because this is included
* by default.)
* @hideinitializer
**********************************************************************/
AZIMUTH = 1U<<9 | CAP_NONE,
/**
* Calculate distance \e s12.
* @hideinitializer
**********************************************************************/
DISTANCE = 1U<<10 | CAP_C1,
/**
* Allow distance \e s12 to be used as input in the direct geodesic
* problem.
* @hideinitializer
**********************************************************************/
DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p,
/**
* Calculate reduced length \e m12.
* @hideinitializer
**********************************************************************/
REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
/**
* Calculate geodesic scales \e M12 and \e M21.
* @hideinitializer
**********************************************************************/
GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
/**
* Calculate area \e S12.
* @hideinitializer
**********************************************************************/
AREA = 1U<<14 | CAP_C4,
/**
* Unroll \e lon2 in the direct calculation.
* @hideinitializer
**********************************************************************/
LONG_UNROLL = 1U<<15,
/**
* All capabilities, calculate everything. (LONG_UNROLL is not
* included in this mask.)
* @hideinitializer
**********************************************************************/
ALL = OUT_ALL| CAP_ALL,
};
/** \name Constructor
**********************************************************************/
///@{
/**
* Constructor for a ellipsoid with
*
* @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
* positive.
**********************************************************************/
Geodesic(real a, real f);
///@}
/** \name Direct geodesic problem specified in terms of distance.
**********************************************************************/
///@{
/**
* Solve the direct geodesic problem where the length of the geodesic
* is specified in terms of distance.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] s12 distance between point 1 and point 2 (meters); it can be
* negative.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
* \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
* 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+. An arc length greater that
* 180&deg; signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shortest
* path: the longitudinal extent must not exceed of 180&deg;.)
*
* The following functions are overloaded versions of Geodesic::Direct
* which omit some of the output parameters. Note, however, that the arc
* length is always computed and returned as the function value.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& m12, real& M12, real& M21, real& S12)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, t, m12, M12, M21, S12);
}
/**
* See the documentation for Geodesic::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE,
lat2, lon2, t, t, t, t, t, t);
}
/**
* See the documentation for Geodesic::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH,
lat2, lon2, azi2, t, t, t, t, t);
}
/**
* See the documentation for Geodesic::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2, real& m12)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
lat2, lon2, azi2, t, m12, t, t, t);
}
/**
* See the documentation for Geodesic::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& M12, real& M21)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
lat2, lon2, azi2, t, t, M12, M21, t);
}
/**
* See the documentation for Geodesic::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& m12, real& M12, real& M21)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, t, m12, M12, M21, t);
}
///@}
/** \name Direct geodesic problem specified in terms of arc length.
**********************************************************************/
///@{
/**
* Solve the direct geodesic problem where the length of the geodesic
* is specified in terms of arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] a12 arc length between point 1 and point 2 (degrees); it can
* be negative.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
* \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
* 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+. An arc length greater that
* 180&deg; signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shortest
* path: the longitudinal extent must not exceed of 180&deg;.)
*
* The following functions are overloaded versions of Geodesic::Direct
* which omit some of the output parameters.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21, real& S12)
const {
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, s12, m12, M12, M21, S12);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE,
lat2, lon2, t, t, t, t, t, t);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH,
lat2, lon2, azi2, t, t, t, t, t);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12)
const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
lat2, lon2, azi2, s12, t, t, t, t);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2,
real& s12, real& m12) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH,
lat2, lon2, azi2, s12, m12, t, t, t);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& M12, real& M21) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
GEODESICSCALE,
lat2, lon2, azi2, s12, t, M12, M21, t);
}
/**
* See the documentation for Geodesic::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, s12, m12, M12, M21, t);
}
///@}
/** \name General version of the direct geodesic solution.
**********************************************************************/
///@{
/**
* The general direct geodesic problem. Geodesic::Direct and
* Geodesic::ArcDirect are defined in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] arcmode boolean flag determining the meaning of the \e
* s12_a12.
* @param[in] s12_a12 if \e arcmode is false, this is the distance between
* point 1 and point 2 (meters); otherwise it is the arc length between
* point 1 and point 2 (degrees); it can be negative.
* @param[in] outmask a bitor'ed combination of Geodesic::mask values
* specifying which of the following parameters should be set.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The Geodesic::mask values possible for \e outmask are
* - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2;
* - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2;
* - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
* - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
* - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
* m12;
* - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21;
* - \e outmask |= Geodesic::AREA for the area \e S12;
* - \e outmask |= Geodesic::ALL for all of the above;
* - \e outmask |= Geodesic::LONG_UNROLL to unroll \e lon2 instead of
* wrapping it into the range [&minus;180&deg;, 180&deg;].
* .
* The function value \e a12 is always computed and returned and this
* equals \e s12_a12 is \e arcmode is true. If \e outmask includes
* Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12.
* It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this
* is automatically included is \e arcmode is false.
*
* With the Geodesic::LONG_UNROLL bit set, the quantity \e lon2 &minus; \e
* lon1 indicates how many times and in what sense the geodesic encircles
* the ellipsoid.
**********************************************************************/
Math::real GenDirect(real lat1, real lon1, real azi1,
bool arcmode, real s12_a12, unsigned outmask,
real& lat2, real& lon2, real& azi2,
real& s12, real& m12, real& M12, real& M21,
real& S12) const;
///@}
/** \name Inverse geodesic problem.
**********************************************************************/
///@{
/**
* Solve the inverse geodesic problem.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
* The values of \e azi1 and \e azi2 returned are in the range
* [&minus;180&deg;, 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+.
*
* The solution to the inverse problem is found using Newton's method. If
* this fails to converge (this is very unlikely in geodetic applications
* but does occur for very eccentric ellipsoids), then the bisection method
* is used to refine the solution.
*
* The following functions are overloaded versions of Geodesic::Inverse
* which omit some of the output parameters. Note, however, that the arc
* length is always computed and returned as the function value.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21, real& S12) const {
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
s12, azi1, azi2, m12, M12, M21, S12);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE,
s12, t, t, t, t, t, t);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& azi1, real& azi2) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
AZIMUTH,
t, azi1, azi2, t, t, t, t);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2)
const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH,
s12, azi1, azi2, t, t, t, t);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12)
const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | REDUCEDLENGTH,
s12, azi1, azi2, m12, t, t, t);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2,
real& M12, real& M21) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | GEODESICSCALE,
s12, azi1, azi2, t, M12, M21, t);
}
/**
* See the documentation for Geodesic::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
s12, azi1, azi2, m12, M12, M21, t);
}
///@}
/** \name General version of inverse geodesic solution.
**********************************************************************/
///@{
/**
* The general inverse geodesic calculation. Geodesic::Inverse is defined
* in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[in] outmask a bitor'ed combination of Geodesic::mask values
* specifying which of the following parameters should be set.
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The Geodesic::mask values possible for \e outmask are
* - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
* - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
* - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
* m12;
* - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21;
* - \e outmask |= Geodesic::AREA for the area \e S12;
* - \e outmask |= Geodesic::ALL for all of the above.
* .
* The arc length is always computed and returned as the function value.
**********************************************************************/
Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
unsigned outmask,
real& s12, real& azi1, real& azi2,
real& m12, real& M12, real& M21, real& S12) const;
///@}
/** \name Interface to GeodesicLine.
**********************************************************************/
///@{
/**
* Set up to compute several points on a single geodesic.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess,
* i.e., which quantities can be returned in calls to
* GeodesicLine::Position.
* @return a GeodesicLine object.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
*
* The Geodesic::mask values are
* - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
* added automatically;
* - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2;
* - \e caps |= Geodesic::AZIMUTH for the azimuth \e azi2; this is
* added automatically;
* - \e caps |= Geodesic::DISTANCE for the distance \e s12;
* - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12;
* - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
* and \e M21;
* - \e caps |= Geodesic::AREA for the area \e S12;
* - \e caps |= Geodesic::DISTANCE_IN permits the length of the
* geodesic to be given in terms of \e s12; without this capability the
* length can only be specified in terms of arc length;
* - \e caps |= Geodesic::ALL for all of the above.
* .
* The default value of \e caps is Geodesic::ALL.
*
* If the point is at a pole, the azimuth is defined by keeping \e lon1
* fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
* limit &epsilon; &rarr; 0+.
**********************************************************************/
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
const;
/**
* Define a GeodesicLine in terms of the inverse geodesic problem.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess,
* i.e., which quantities can be returned in calls to
* GeodesicLine::Position.
* @return a GeodesicLine object.
*
* This function sets point 3 of the GeodesicLine to correspond to point 2
* of the inverse geodesic problem.
*
* \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLine InverseLine(real lat1, real lon1, real lat2, real lon2,
unsigned caps = ALL) const;
/**
* Define a GeodesicLine in terms of the direct geodesic problem specified
* in terms of distance.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] s12 distance between point 1 and point 2 (meters); it can be
* negative.
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess,
* i.e., which quantities can be returned in calls to
* GeodesicLine::Position.
* @return a GeodesicLine object.
*
* This function sets point 3 of the GeodesicLine to correspond to point 2
* of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLine DirectLine(real lat1, real lon1, real azi1, real s12,
unsigned caps = ALL) const;
/**
* Define a GeodesicLine in terms of the direct geodesic problem specified
* in terms of arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] a12 arc length between point 1 and point 2 (degrees); it can
* be negative.
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess,
* i.e., which quantities can be returned in calls to
* GeodesicLine::Position.
* @return a GeodesicLine object.
*
* This function sets point 3 of the GeodesicLine to correspond to point 2
* of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLine ArcDirectLine(real lat1, real lon1, real azi1, real a12,
unsigned caps = ALL) const;
/**
* Define a GeodesicLine in terms of the direct geodesic problem specified
* in terms of either distance or arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] arcmode boolean flag determining the meaning of the \e
* s12_a12.
* @param[in] s12_a12 if \e arcmode is false, this is the distance between
* point 1 and point 2 (meters); otherwise it is the arc length between
* point 1 and point 2 (degrees); it can be negative.
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess,
* i.e., which quantities can be returned in calls to
* GeodesicLine::Position.
* @return a GeodesicLine object.
*
* This function sets point 3 of the GeodesicLine to correspond to point 2
* of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLine GenDirectLine(real lat1, real lon1, real azi1,
bool arcmode, real s12_a12,
unsigned caps = ALL) const;
///@}
/** \name Inspector functions.
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _a; }
/**
* @return \e f the flattening of the ellipsoid. This is the
* value used in the constructor.
**********************************************************************/
Math::real Flattening() const { return _f; }
/**
* @return total area of ellipsoid in meters<sup>2</sup>. The area of a
* polygon encircling a pole can be found by adding
* Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
* polygon.
**********************************************************************/
Math::real EllipsoidArea() const
{ return 4 * Math::pi() * _c2; }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
/**
* A global instantiation of Geodesic with the parameters for the WGS84
* ellipsoid.
**********************************************************************/
static const Geodesic& WGS84();
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_GEODESIC_HPP

View File

@@ -0,0 +1,869 @@
/**
* \file GeodesicExact.hpp
* \brief Header for GeographicLib::GeodesicExact class
*
* Copyright (c) Charles Karney (2012-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESICEXACT_HPP)
#define GEOGRAPHICLIB_GEODESICEXACT_HPP 1
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/EllipticFunction.hpp>
#if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER)
/**
* The order of the expansions used by GeodesicExact.
**********************************************************************/
# define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30
#endif
namespace GeographicLib {
class GeodesicLineExact;
/**
* \brief Exact geodesic calculations
*
* The equations for geodesics on an ellipsoid can be expressed in terms of
* incomplete elliptic integrals. The Geodesic class expands these integrals
* in a series in the flattening \e f and this provides an accurate solution
* for \e f &isin; [-0.01, 0.01]. The GeodesicExact class computes the
* ellitpic integrals directly and so provides a solution which is valid for
* all \e f. However, in practice, its use should be limited to about
* <i>b</i>/\e a &isin; [0.01, 100] or \e f &isin; [&minus;99, 0.99].
*
* For the WGS84 ellipsoid, these classes are 2--3 times \e slower than the
* series solution and 2--3 times \e less \e accurate (because it's less easy
* to control round-off errors with the elliptic integral formulation); i.e.,
* the error is about 40 nm (40 nanometers) instead of 15 nm. However the
* error in the series solution scales as <i>f</i><sup>7</sup> while the
* error in the elliptic integral solution depends weakly on \e f. If the
* quarter meridian distance is 10000 km and the ratio <i>b</i>/\e a = 1
* &minus; \e f is varied then the approximate maximum error (expressed as a
* distance) is <pre>
* 1 - f error (nm)
* 1/128 387
* 1/64 345
* 1/32 269
* 1/16 210
* 1/8 115
* 1/4 69
* 1/2 36
* 1 15
* 2 25
* 4 96
* 8 318
* 16 985
* 32 2352
* 64 6008
* 128 19024
* </pre>
*
* The computation of the area in these classes is via a 30th order series.
* This gives accurate results for <i>b</i>/\e a &isin; [1/2, 2]; the
* accuracy is about 8 decimal digits for <i>b</i>/\e a &isin; [1/4, 4].
*
* See \ref geodellip for the formulation. See the documentation on the
* Geodesic class for additional information on the geodesic problems.
*
* Example of use:
* \include example-GeodesicExact.cpp
*
* <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
* providing access to the functionality of GeodesicExact and
* GeodesicLineExact (via the -E option).
**********************************************************************/
class GEOGRAPHICLIB_EXPORT GeodesicExact {
private:
typedef Math::real real;
friend class GeodesicLineExact;
static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
static const unsigned maxit1_ = 20;
unsigned maxit2_;
real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
enum captype {
CAP_NONE = 0U,
CAP_E = 1U<<0,
// Skip 1U<<1 for compatibility with Geodesic (not required)
CAP_D = 1U<<2,
CAP_H = 1U<<3,
CAP_C4 = 1U<<4,
CAP_ALL = 0x1FU,
CAP_MASK = CAP_ALL,
OUT_ALL = 0x7F80U,
OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
};
static real CosSeries(real sinx, real cosx, const real c[], int n);
static real Astroid(real x, real y);
real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
real _C4x[nC4x_];
void Lengths(const EllipticFunction& E,
real sig12,
real ssig1, real csig1, real dn1,
real ssig2, real csig2, real dn2,
real cbet1, real cbet2, unsigned outmask,
real& s12s, real& m12a, real& m0,
real& M12, real& M21) const;
real InverseStart(EllipticFunction& E,
real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real lam12, real slam12, real clam12,
real& salp1, real& calp1,
real& salp2, real& calp2, real& dnm) const;
real Lambda12(real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real salp1, real calp1, real slam120, real clam120,
real& salp2, real& calp2, real& sig12,
real& ssig1, real& csig1, real& ssig2, real& csig2,
EllipticFunction& E,
real& domg12, bool diffp, real& dlam12) const;
real GenInverse(real lat1, real lon1, real lat2, real lon2,
unsigned outmask, real& s12,
real& salp1, real& calp1, real& salp2, real& calp2,
real& m12, real& M12, real& M21, real& S12) const;
// These are Maxima generated functions to provide series approximations to
// the integrals for the area.
void C4coeff();
void C4f(real k2, real c[]) const;
// Large coefficients are split so that lo contains the low 52 bits and hi
// the rest. This choice avoids double rounding with doubles and higher
// precision types. float coefficients will suffer double rounding;
// however the accuracy is already lousy for floats.
static Math::real reale(long long hi, long long lo) {
using std::ldexp;
return ldexp(real(hi), 52) + lo;
}
public:
/**
* Bit masks for what calculations to do. These masks do double duty.
* They signify to the GeodesicLineExact::GeodesicLineExact constructor and
* to GeodesicExact::Line what capabilities should be included in the
* GeodesicLineExact object. They also specify which results to return in
* the general routines GeodesicExact::GenDirect and
* GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a
* duplication of this enum.
**********************************************************************/
enum mask {
/**
* No capabilities, no output.
* @hideinitializer
**********************************************************************/
NONE = 0U,
/**
* Calculate latitude \e lat2. (It's not necessary to include this as a
* capability to GeodesicLineExact because this is included by default.)
* @hideinitializer
**********************************************************************/
LATITUDE = 1U<<7 | CAP_NONE,
/**
* Calculate longitude \e lon2.
* @hideinitializer
**********************************************************************/
LONGITUDE = 1U<<8 | CAP_H,
/**
* Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
* include this as a capability to GeodesicLineExact because this is
* included by default.)
* @hideinitializer
**********************************************************************/
AZIMUTH = 1U<<9 | CAP_NONE,
/**
* Calculate distance \e s12.
* @hideinitializer
**********************************************************************/
DISTANCE = 1U<<10 | CAP_E,
/**
* Allow distance \e s12 to be used as input in the direct geodesic
* problem.
* @hideinitializer
**********************************************************************/
DISTANCE_IN = 1U<<11 | CAP_E,
/**
* Calculate reduced length \e m12.
* @hideinitializer
**********************************************************************/
REDUCEDLENGTH = 1U<<12 | CAP_D,
/**
* Calculate geodesic scales \e M12 and \e M21.
* @hideinitializer
**********************************************************************/
GEODESICSCALE = 1U<<13 | CAP_D,
/**
* Calculate area \e S12.
* @hideinitializer
**********************************************************************/
AREA = 1U<<14 | CAP_C4,
/**
* Unroll \e lon2 in the direct calculation.
* @hideinitializer
**********************************************************************/
LONG_UNROLL = 1U<<15,
/**
* All capabilities, calculate everything. (LONG_UNROLL is not
* included in this mask.)
* @hideinitializer
**********************************************************************/
ALL = OUT_ALL| CAP_ALL,
};
/** \name Constructor
**********************************************************************/
///@{
/**
* Constructor for a ellipsoid with
*
* @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
* Negative \e f gives a prolate ellipsoid.
* @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
* positive.
**********************************************************************/
GeodesicExact(real a, real f);
///@}
/** \name Direct geodesic problem specified in terms of distance.
**********************************************************************/
///@{
/**
* Perform the direct geodesic calculation where the length of the geodesic
* is specified in terms of distance.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] s12 distance between point 1 and point 2 (meters); it can be
* signed.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
* \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
* 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+. An arc length greater that
* 180&deg; signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shortest
* path: the longitudinal extent must not exceed of 180&deg;.)
*
* The following functions are overloaded versions of GeodesicExact::Direct
* which omit some of the output parameters. Note, however, that the arc
* length is always computed and returned as the function value.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& m12, real& M12, real& M21, real& S12)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, t, m12, M12, M21, S12);
}
/**
* See the documentation for GeodesicExact::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE,
lat2, lon2, t, t, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH,
lat2, lon2, azi2, t, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2, real& m12)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
lat2, lon2, azi2, t, m12, t, t, t);
}
/**
* See the documentation for GeodesicExact::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& M12, real& M21)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
lat2, lon2, azi2, t, t, M12, M21, t);
}
/**
* See the documentation for GeodesicExact::Direct.
**********************************************************************/
Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& m12, real& M12, real& M21)
const {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, t, m12, M12, M21, t);
}
///@}
/** \name Direct geodesic problem specified in terms of arc length.
**********************************************************************/
///@{
/**
* Perform the direct geodesic calculation where the length of the geodesic
* is specified in terms of arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] a12 arc length between point 1 and point 2 (degrees); it can
* be signed.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
* \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
* 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+. An arc length greater that
* 180&deg; signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shortest
* path: the longitudinal extent must not exceed of 180&deg;.)
*
* The following functions are overloaded versions of GeodesicExact::Direct
* which omit some of the output parameters.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21, real& S12)
const {
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, s12, m12, M12, M21, S12);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE,
lat2, lon2, t, t, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH,
lat2, lon2, azi2, t, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12)
const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
lat2, lon2, azi2, s12, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2,
real& s12, real& m12) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH,
lat2, lon2, azi2, s12, m12, t, t, t);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& M12, real& M21) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
GEODESICSCALE,
lat2, lon2, azi2, s12, t, M12, M21, t);
}
/**
* See the documentation for GeodesicExact::ArcDirect.
**********************************************************************/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21) const {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, s12, m12, M12, M21, t);
}
///@}
/** \name General version of the direct geodesic solution.
**********************************************************************/
///@{
/**
* The general direct geodesic calculation. GeodesicExact::Direct and
* GeodesicExact::ArcDirect are defined in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] arcmode boolean flag determining the meaning of the second
* parameter.
* @param[in] s12_a12 if \e arcmode is false, this is the distance between
* point 1 and point 2 (meters); otherwise it is the arc length between
* point 1 and point 2 (degrees); it can be signed.
* @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
* specifying which of the following parameters should be set.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The GeodesicExact::mask values possible for \e outmask are
* - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
* - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
* - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
* - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
* - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
* m12;
* - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21;
* - \e outmask |= GeodesicExact::AREA for the area \e S12;
* - \e outmask |= GeodesicExact::ALL for all of the above;
* - \e outmask |= GeodesicExact::LONG_UNROLL to unroll \e lon2 instead of
* wrapping it into the range [&minus;180&deg;, 180&deg;].
* .
* The function value \e a12 is always computed and returned and this
* equals \e s12_a12 is \e arcmode is true. If \e outmask includes
* GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
* s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
* \e outmask; this is automatically included is \e arcmode is false.
*
* With the GeodesicExact::LONG_UNROLL bit set, the quantity \e lon2
* &minus; \e lon1 indicates how many times and in what sense the geodesic
* encircles the ellipsoid.
**********************************************************************/
Math::real GenDirect(real lat1, real lon1, real azi1,
bool arcmode, real s12_a12, unsigned outmask,
real& lat2, real& lon2, real& azi2,
real& s12, real& m12, real& M12, real& M21,
real& S12) const;
///@}
/** \name Inverse geodesic problem.
**********************************************************************/
///@{
/**
* Perform the inverse geodesic calculation.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
* The values of \e azi1 and \e azi2 returned are in the range
* [&minus;180&deg;, 180&deg;].
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
* and taking the limit &epsilon; &rarr; 0+.
*
* The following functions are overloaded versions of
* GeodesicExact::Inverse which omit some of the output parameters. Note,
* however, that the arc length is always computed and returned as the
* function value.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21, real& S12) const {
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
s12, azi1, azi2, m12, M12, M21, S12);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE,
s12, t, t, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& azi1, real& azi2) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
AZIMUTH,
t, azi1, azi2, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2)
const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH,
s12, azi1, azi2, t, t, t, t);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12)
const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | REDUCEDLENGTH,
s12, azi1, azi2, m12, t, t, t);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2,
real& M12, real& M21) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | GEODESICSCALE,
s12, azi1, azi2, t, M12, M21, t);
}
/**
* See the documentation for GeodesicExact::Inverse.
**********************************************************************/
Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21) const {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
s12, azi1, azi2, m12, M12, M21, t);
}
///@}
/** \name General version of inverse geodesic solution.
**********************************************************************/
///@{
/**
* The general inverse geodesic calculation. GeodesicExact::Inverse is
* defined in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
* specifying which of the following parameters should be set.
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The GeodesicExact::mask values possible for \e outmask are
* - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
* - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
* - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
* m12;
* - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21;
* - \e outmask |= GeodesicExact::AREA for the area \e S12;
* - \e outmask |= GeodesicExact::ALL for all of the above.
* .
* The arc length is always computed and returned as the function value.
**********************************************************************/
Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
unsigned outmask,
real& s12, real& azi1, real& azi2,
real& m12, real& M12, real& M21, real& S12) const;
///@}
/** \name Interface to GeodesicLineExact.
**********************************************************************/
///@{
/**
* Set up to compute several points on a single geodesic.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] caps bitor'ed combination of GeodesicExact::mask values
* specifying the capabilities the GeodesicLineExact object should
* possess, i.e., which quantities can be returned in calls to
* GeodesicLineExact::Position.
* @return a GeodesicLineExact object.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
*
* The GeodesicExact::mask values are
* - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is
* added automatically;
* - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2;
* - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is
* added automatically;
* - \e caps |= GeodesicExact::DISTANCE for the distance \e s12;
* - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12;
* - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12
* and \e M21;
* - \e caps |= GeodesicExact::AREA for the area \e S12;
* - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the
* geodesic to be given in terms of \e s12; without this capability the
* length can only be specified in terms of arc length;
* - \e caps |= GeodesicExact::ALL for all of the above.
* .
* The default value of \e caps is GeodesicExact::ALL which turns on all
* the capabilities.
*
* If the point is at a pole, the azimuth is defined by keeping \e lon1
* fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
* limit &epsilon; &rarr; 0+.
**********************************************************************/
GeodesicLineExact Line(real lat1, real lon1, real azi1,
unsigned caps = ALL) const;
/**
* Define a GeodesicLineExact in terms of the inverse geodesic problem.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[in] caps bitor'ed combination of GeodesicExact::mask values
* specifying the capabilities the GeodesicLineExact object should
* possess, i.e., which quantities can be returned in calls to
* GeodesicLineExact::Position.
* @return a GeodesicLineExact object.
*
* This function sets point 3 of the GeodesicLineExact to correspond to
* point 2 of the inverse geodesic problem.
*
* \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2,
unsigned caps = ALL) const;
/**
* Define a GeodesicLineExact in terms of the direct geodesic problem
* specified in terms of distance.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] s12 distance between point 1 and point 2 (meters); it can be
* negative.
* @param[in] caps bitor'ed combination of GeodesicExact::mask values
* specifying the capabilities the GeodesicLineExact object should
* possess, i.e., which quantities can be returned in calls to
* GeodesicLineExact::Position.
* @return a GeodesicLineExact object.
*
* This function sets point 3 of the GeodesicLineExact to correspond to
* point 2 of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12,
unsigned caps = ALL) const;
/**
* Define a GeodesicLineExact in terms of the direct geodesic problem
* specified in terms of arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] a12 arc length between point 1 and point 2 (degrees); it can
* be negative.
* @param[in] caps bitor'ed combination of GeodesicExact::mask values
* specifying the capabilities the GeodesicLineExact object should
* possess, i.e., which quantities can be returned in calls to
* GeodesicLineExact::Position.
* @return a GeodesicLineExact object.
*
* This function sets point 3 of the GeodesicLineExact to correspond to
* point 2 of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12,
unsigned caps = ALL) const;
/**
* Define a GeodesicLineExact in terms of the direct geodesic problem
* specified in terms of either distance or arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] arcmode boolean flag determining the meaning of the \e
* s12_a12.
* @param[in] s12_a12 if \e arcmode is false, this is the distance between
* point 1 and point 2 (meters); otherwise it is the arc length between
* point 1 and point 2 (degrees); it can be negative.
* @param[in] caps bitor'ed combination of GeodesicExact::mask values
* specifying the capabilities the GeodesicLineExact object should
* possess, i.e., which quantities can be returned in calls to
* GeodesicLineExact::Position.
* @return a GeodesicLineExact object.
*
* This function sets point 3 of the GeodesicLineExact to correspond to
* point 2 of the direct geodesic problem.
*
* \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
**********************************************************************/
GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1,
bool arcmode, real s12_a12,
unsigned caps = ALL) const;
///@}
/** \name Inspector functions.
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor.
**********************************************************************/
Math::real EquatorialRadius() const { return _a; }
/**
* @return \e f the flattening of the ellipsoid. This is the
* value used in the constructor.
**********************************************************************/
Math::real Flattening() const { return _f; }
/**
* @return total area of ellipsoid in meters<sup>2</sup>. The area of a
* polygon encircling a pole can be found by adding
* GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
* the polygon.
**********************************************************************/
Math::real EllipsoidArea() const
{ return 4 * Math::pi() * _c2; }
/**
* \deprecated An old name for EquatorialRadius().
**********************************************************************/
GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
Math::real MajorRadius() const { return EquatorialRadius(); }
///@}
/**
* A global instantiation of GeodesicExact with the parameters for the
* WGS84 ellipsoid.
**********************************************************************/
static const GeodesicExact& WGS84();
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_GEODESICEXACT_HPP

View File

@@ -0,0 +1,708 @@
/**
* \file GeodesicLine.hpp
* \brief Header for GeographicLib::GeodesicLine class
*
* Copyright (c) Charles Karney (2009-2020) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP)
#define GEOGRAPHICLIB_GEODESICLINE_HPP 1
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Geodesic.hpp>
namespace GeographicLib {
/**
* \brief A geodesic line
*
* GeodesicLine facilitates the determination of a series of points on a
* single geodesic. The starting point (\e lat1, \e lon1) and the azimuth \e
* azi1 are specified in the constructor; alternatively, the Geodesic::Line
* method can be used to create a GeodesicLine. GeodesicLine.Position
* returns the location of point 2 a distance \e s12 along the geodesic. In
* addition, GeodesicLine.ArcPosition gives the position of point 2 an arc
* length \e a12 along the geodesic.
*
* You can register the position of a reference point 3 a distance (arc
* length), \e s13 (\e a13) along the geodesic with the
* GeodesicLine.SetDistance (GeodesicLine.SetArc) functions. Points a
* fractional distance along the line can be found by providing, for example,
* 0.5 * Distance() as an argument to GeodesicLine.Position. The
* Geodesic::InverseLine or Geodesic::DirectLine methods return GeodesicLine
* objects with point 3 set to the point 2 of the corresponding geodesic
* problem. GeodesicLine objects created with the public constructor or with
* Geodesic::Line have \e s13 and \e a13 set to NaNs.
*
* The default copy constructor and assignment operators work with this
* class. Similarly, a vector can be used to hold GeodesicLine objects.
*
* The calculations are accurate to better than 15 nm (15 nanometers). See
* Sec. 9 of
* <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for
* details. The algorithms used by this class are based on series expansions
* using the flattening \e f as a small parameter. These are only accurate
* for |<i>f</i>| &lt; 0.02; however reasonably accurate results will be
* obtained for |<i>f</i>| &lt; 0.2. For very eccentric ellipsoids, use
* GeodesicLineExact instead.
*
* The algorithms are described in
* - C. F. F. Karney,
* <a href="https://doi.org/10.1007/s00190-012-0578-z">
* Algorithms for geodesics</a>,
* J. Geodesy <b>87</b>, 43--55 (2013);
* DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
* 10.1007/s00190-012-0578-z</a>;
* addenda:
* <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
* geod-addenda.html</a>.
* .
* For more information on geodesics see \ref geodesic.
*
* Example of use:
* \include example-GeodesicLine.cpp
*
* <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
* providing access to the functionality of Geodesic and GeodesicLine.
**********************************************************************/
class GEOGRAPHICLIB_EXPORT GeodesicLine {
private:
typedef Math::real real;
friend class Geodesic;
static const int nC1_ = Geodesic::nC1_;
static const int nC1p_ = Geodesic::nC1p_;
static const int nC2_ = Geodesic::nC2_;
static const int nC3_ = Geodesic::nC3_;
static const int nC4_ = Geodesic::nC4_;
real tiny_;
real _lat1, _lon1, _azi1;
real _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
_salp1, _calp1, _ssig1, _csig1, _dn1, _stau1, _ctau1, _somg1, _comg1,
_A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41;
real _a13, _s13;
// index zero elements of _C1a, _C1pa, _C2a, _C3a are unused
real _C1a[nC1_ + 1], _C1pa[nC1p_ + 1], _C2a[nC2_ + 1], _C3a[nC3_],
_C4a[nC4_]; // all the elements of _C4a are used
unsigned _caps;
void LineInit(const Geodesic& g,
real lat1, real lon1,
real azi1, real salp1, real calp1,
unsigned caps);
GeodesicLine(const Geodesic& g,
real lat1, real lon1,
real azi1, real salp1, real calp1,
unsigned caps, bool arcmode, real s13_a13);
enum captype {
CAP_NONE = Geodesic::CAP_NONE,
CAP_C1 = Geodesic::CAP_C1,
CAP_C1p = Geodesic::CAP_C1p,
CAP_C2 = Geodesic::CAP_C2,
CAP_C3 = Geodesic::CAP_C3,
CAP_C4 = Geodesic::CAP_C4,
CAP_ALL = Geodesic::CAP_ALL,
CAP_MASK = Geodesic::CAP_MASK,
OUT_ALL = Geodesic::OUT_ALL,
OUT_MASK = Geodesic::OUT_MASK,
};
public:
/**
* Bit masks for what calculations to do. They signify to the
* GeodesicLine::GeodesicLine constructor and to Geodesic::Line what
* capabilities should be included in the GeodesicLine object. This is
* merely a duplication of Geodesic::mask.
**********************************************************************/
enum mask {
/**
* No capabilities, no output.
* @hideinitializer
**********************************************************************/
NONE = Geodesic::NONE,
/**
* Calculate latitude \e lat2. (It's not necessary to include this as a
* capability to GeodesicLine because this is included by default.)
* @hideinitializer
**********************************************************************/
LATITUDE = Geodesic::LATITUDE,
/**
* Calculate longitude \e lon2.
* @hideinitializer
**********************************************************************/
LONGITUDE = Geodesic::LONGITUDE,
/**
* Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
* include this as a capability to GeodesicLine because this is included
* by default.)