1 // Boost.Geometry - gis-projections (based on PROJ4) 2 3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. 4 5 // This file was modified by Oracle on 2017, 2018, 2019. 6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates. 7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. 8 9 // Use, modification and distribution is subject to the Boost Software License, 10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 11 // http://www.boost.org/LICENSE_1_0.txt) 12 13 // This file is converted from PROJ4, http://trac.osgeo.org/proj 14 // PROJ4 is originally written by Gerald Evenden (then of the USGS) 15 // PROJ4 is maintained by Frank Warmerdam 16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels 17 18 // Last updated version of proj: 5.0.0 19 20 // Original copyright notice: 21 22 // Permission is hereby granted, free of charge, to any person obtaining a 23 // copy of this software and associated documentation files (the "Software"), 24 // to deal in the Software without restriction, including without limitation 25 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 26 // and/or sell copies of the Software, and to permit persons to whom the 27 // Software is furnished to do so, subject to the following conditions: 28 29 // The above copyright notice and this permission notice shall be included 30 // in all copies or substantial portions of the Software. 31 32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 38 // DEALINGS IN THE SOFTWARE. 39 40 #ifndef BOOST_GEOMETRY_PROJECTIONS_CASS_HPP 41 #define BOOST_GEOMETRY_PROJECTIONS_CASS_HPP 42 43 #include <boost/geometry/srs/projections/impl/base_static.hpp> 44 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 45 #include <boost/geometry/srs/projections/impl/projects.hpp> 46 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 47 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp> 48 49 50 namespace boost { namespace geometry 51 { 52 53 namespace projections 54 { 55 #ifndef DOXYGEN_NO_DETAIL 56 namespace detail { namespace cass 57 { 58 59 //static const double epsilon10 = 1e-10; 60 //static const double C1 = .16666666666666666666; 61 //static const double C2 = .00833333333333333333; 62 //static const double C3 = .04166666666666666666; 63 //static const double C4 = .33333333333333333333; 64 //static const double C5 = .06666666666666666666; 65 66 template <typename T> C1()67 inline T C1() { return .16666666666666666666666666666666666666; } 68 template <typename T> C2()69 inline T C2() { return .00833333333333333333333333333333333333; } 70 template <typename T> C3()71 inline T C3() { return .04166666666666666666666666666666666666; } 72 template <typename T> C4()73 inline T C4() { return .33333333333333333333333333333333333333; } 74 template <typename T> C5()75 inline T C5() { return .06666666666666666666666666666666666666; } 76 77 template <typename T> 78 struct par_cass 79 { 80 T m0; 81 detail::en<T> en; 82 }; 83 84 template <typename T, typename Parameters> 85 struct base_cass_ellipsoid 86 { 87 par_cass<T> m_proj_parm; 88 89 // FORWARD(e_forward) ellipsoid 90 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::cass::base_cass_ellipsoid91 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 92 { 93 static const T C1 = cass::C1<T>(); 94 static const T C2 = cass::C2<T>(); 95 static const T C3 = cass::C3<T>(); 96 97 T n = sin(lp_lat); 98 T c = cos(lp_lat); 99 xy_y = pj_mlfn(lp_lat, n, c, this->m_proj_parm.en); 100 n = 1./sqrt(1. - par.es * n * n); 101 T tn = tan(lp_lat); 102 T t = tn * tn; 103 T a1 = lp_lon * c; 104 c *= par.es * c / (1 - par.es); 105 T a2 = a1 * a1; 106 xy_x = n * a1 * (1. - a2 * t * 107 (C1 - (8. - t + 8. * c) * a2 * C2)); 108 xy_y -= this->m_proj_parm.m0 - n * tn * a2 * 109 (.5 + (5. - t + 6. * c) * a2 * C3); 110 } 111 112 // INVERSE(e_inverse) ellipsoid 113 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::cass::base_cass_ellipsoid114 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 115 { 116 static const T C3 = cass::C3<T>(); 117 static const T C4 = cass::C4<T>(); 118 static const T C5 = cass::C5<T>(); 119 120 T ph1; 121 122 ph1 = pj_inv_mlfn(this->m_proj_parm.m0 + xy_y, par.es, this->m_proj_parm.en); 123 T tn = tan(ph1); T t = tn * tn; 124 T n = sin(ph1); 125 T r = 1. / (1. - par.es * n * n); 126 n = sqrt(r); 127 r *= (1. - par.es) * n; 128 T dd = xy_x / n; 129 T d2 = dd * dd; 130 lp_lat = ph1 - (n * tn / r) * d2 * 131 (.5 - (1. + 3. * t) * d2 * C3); 132 lp_lon = dd * (1. + t * d2 * 133 (-C4 + (1. + 3. * t) * d2 * C5)) / cos(ph1); 134 } 135 get_nameboost::geometry::projections::detail::cass::base_cass_ellipsoid136 static inline std::string get_name() 137 { 138 return "cass_ellipsoid"; 139 } 140 141 }; 142 143 template <typename T, typename Parameters> 144 struct base_cass_spheroid 145 { 146 par_cass<T> m_proj_parm; 147 148 // FORWARD(s_forward) spheroid 149 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::cass::base_cass_spheroid150 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 151 { 152 xy_x = asin(cos(lp_lat) * sin(lp_lon)); 153 xy_y = atan2(tan(lp_lat) , cos(lp_lon)) - par.phi0; 154 } 155 156 // INVERSE(s_inverse) spheroid 157 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::cass::base_cass_spheroid158 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 159 { 160 T dd = xy_y + par.phi0; 161 lp_lat = asin(sin(dd) * cos(xy_x)); 162 lp_lon = atan2(tan(xy_x), cos(dd)); 163 } 164 get_nameboost::geometry::projections::detail::cass::base_cass_spheroid165 static inline std::string get_name() 166 { 167 return "cass_spheroid"; 168 } 169 170 }; 171 172 // Cassini 173 template <typename Parameters, typename T> setup_cass(Parameters & par,par_cass<T> & proj_parm)174 inline void setup_cass(Parameters& par, par_cass<T>& proj_parm) 175 { 176 if (par.es) { 177 proj_parm.en = pj_enfn<T>(par.es); 178 proj_parm.m0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en); 179 } else { 180 } 181 } 182 183 }} // namespace detail::cass 184 #endif // doxygen 185 186 /*! 187 \brief Cassini projection 188 \ingroup projections 189 \tparam Geographic latlong point type 190 \tparam Cartesian xy point type 191 \tparam Parameters parameter type 192 \par Projection characteristics 193 - Cylindrical 194 - Spheroid 195 - Ellipsoid 196 \par Example 197 \image html ex_cass.gif 198 */ 199 template <typename T, typename Parameters> 200 struct cass_ellipsoid : public detail::cass::base_cass_ellipsoid<T, Parameters> 201 { 202 template <typename Params> cass_ellipsoidboost::geometry::projections::cass_ellipsoid203 inline cass_ellipsoid(Params const& , Parameters const& par) 204 { 205 detail::cass::setup_cass(par, this->m_proj_parm); 206 } 207 }; 208 209 /*! 210 \brief Cassini projection 211 \ingroup projections 212 \tparam Geographic latlong point type 213 \tparam Cartesian xy point type 214 \tparam Parameters parameter type 215 \par Projection characteristics 216 - Cylindrical 217 - Spheroid 218 - Ellipsoid 219 \par Example 220 \image html ex_cass.gif 221 */ 222 template <typename T, typename Parameters> 223 struct cass_spheroid : public detail::cass::base_cass_spheroid<T, Parameters> 224 { 225 template <typename Params> cass_spheroidboost::geometry::projections::cass_spheroid226 inline cass_spheroid(Params const& , Parameters const& par) 227 { 228 detail::cass::setup_cass(par, this->m_proj_parm); 229 } 230 }; 231 232 #ifndef DOXYGEN_NO_DETAIL 233 namespace detail 234 { 235 236 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_cass,cass_spheroid,cass_ellipsoid)237 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_cass, cass_spheroid, cass_ellipsoid) 238 239 // Factory entry(s) 240 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(cass_entry, cass_spheroid, cass_ellipsoid) 241 242 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(cass_init) 243 { 244 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(cass, cass_entry); 245 } 246 247 } // namespace detail 248 #endif // doxygen 249 250 } // namespace projections 251 252 }} // namespace boost::geometry 253 254 #endif // BOOST_GEOMETRY_PROJECTIONS_CASS_HPP 255 256