File indexing completed on 2024-04-21 14:50:55

0001 // SPDX-License-Identifier: LGPL-2.1-or-later
0002 //
0003 // SPDX-FileCopyrightText: 2007 Torsten Rahn <tackat@kde.org>
0004 // SPDX-FileCopyrightText: 2007 Inge Wallin <ingwa@kde.org>
0005 //
0006 
0007 #ifndef MARBLE_MARBLEMATH_H
0008 #define MARBLE_MARBLEMATH_H
0009 
0010 #include <QtGlobal>
0011 
0012 #include <cmath>
0013 
0014 
0015 namespace
0016 {
0017     const qreal a1 = 1.0/6.0; 
0018     const qreal a2 = 1.0/24.0; 
0019     const qreal a3 = 61.0/5040; 
0020     const qreal a4 = 277.0/72576.0;  
0021     const qreal a5 = 50521.0/39916800.0; 
0022     const qreal a6 = 41581.0/95800320.0; 
0023     const qreal a7 = 199360981.0/1307674368000.0; 
0024     const qreal a8 = 228135437.0/4184557977600.0; 
0025     const qreal a9 = 2404879675441.0/121645100408832000.0; 
0026     const qreal a10 = 14814847529501.0/2043637686868377600.0; 
0027     const qreal a11 = 69348874393137901.0/25852016738884976640000.0; 
0028     const qreal a12 = 238685140977801337.0/238634000666630553600000.0; 
0029     const qreal a13 = 4087072509293123892361.0/10888869450418352160768000000.0;
0030     const qreal a14 = 454540704683713199807.0/3209350995912777478963200000.0;
0031     const qreal a15 = 441543893249023104553682821.0/8222838654177922817725562880000000.0;
0032     const qreal a16 = 2088463430347521052196056349.0/102156677868375135241390522368000000.0;
0033 }
0034 
0035 namespace Marble
0036 {
0037 
0038 /**
0039  * @brief This method calculates the shortest distance between two points on a sphere.
0040  * @brief See: https://en.wikipedia.org/wiki/Great-circle_distance
0041  * @param lon1 longitude of first point in radians
0042  * @param lat1 latitude of first point in radians
0043  * @param lon2 longitude of second point in radians
0044  * @param lat2 latitude of second point in radians
0045  */
0046 inline qreal distanceSphere( qreal lon1, qreal lat1, qreal lon2, qreal lat2 ) {
0047 
0048     qreal h1 = sin( 0.5 * ( lat2 - lat1 ) );
0049     qreal h2 = sin( 0.5 * ( lon2 - lon1 ) );
0050     qreal d = h1 * h1 + cos( lat1 ) * cos( lat2 ) * h2 * h2;
0051 
0052     return 2.0 * atan2( sqrt( d ), sqrt( 1.0 - d ) );
0053 }
0054 
0055 
0056 /**
0057  * @brief This method roughly calculates the shortest distance between two points on a sphere.
0058  * @brief It's probably faster than distanceSphere(...) but for 7 significant digits only has
0059  * @brief an accuracy of about 1 arcmin.
0060  * @brief See: https://en.wikipedia.org/wiki/Great-circle_distance
0061  */
0062 inline qreal distanceSphereApprox( qreal lon1, qreal lat1, qreal lon2, qreal lat2 ) {
0063     return acos( sin( lat1 ) * sin( lat2 ) + cos( lat1 ) * cos( lat2 ) * cos( lon1 - lon2 ) );
0064 }
0065 
0066 
0067 /**
0068  * @brief This method is a fast Mac Laurin power series approximation of the 
0069  * @brief inverse Gudermannian. The inverse Gudermannian gives the vertical 
0070  * @brief position y in the Mercator projection in terms of the latitude.
0071  * @brief See: https://en.wikipedia.org/wiki/Mercator_projection
0072  */
0073 inline qreal gdInv( qreal x ) {
0074         const qreal x2 = x * x;
0075         return x 
0076             + x * x2 * (  a1
0077             + x2 * ( a2  + x2 * ( a3  + x2 * ( a4  + x2 * ( a5
0078             + x2 * ( a6  + x2 * ( a7  + x2 * ( a8  + x2 * ( a9
0079             + x2 * ( a10 + x2 * ( a11 + x2 * ( a12 + x2 * ( a13
0080             + x2 * ( a14 + x2 * ( a15 + x2 * ( a16 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) );
0081 }
0082 
0083 }
0084 
0085 /**
0086  * @brief This method is a fast Mac Laurin power series approximation of the
0087  *        Gudermannian. The Gudermannian gives the latitude
0088  *        in the Mercator projection in terms of the vertical position y.
0089  *        See: https://en.wikipedia.org/wiki/Mercator_projection
0090  */
0091 inline qreal gd( qreal x ) {
0092 
0093     /*
0094     const qreal x2 = x * x;
0095     return x
0096          - x * x2 * (  a1
0097          - x2 * ( a2  - x2 * ( a3  - x2 * ( a4  - x2 * ( a5
0098          - x2 * ( a6  - x2 * ( a7  - x2 * ( a8  - x2 * ( a9
0099          - x2 * ( a10 - x2 * ( a11 - x2 * ( a12 - x2 * ( a13
0100          - x2 * ( a14 - x2 * ( a15 - x2 * ( a16 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) );
0101     */
0102 
0103     return atan ( sinh ( x ) );
0104 }
0105 
0106 #endif