← Back to branch summary

~ubuntu-branches/ubuntu/saucy/merkaartor/saucy

« back to all changes in this revision

Viewing changes to include/builtin-ggl/ggl/extensions/gis/projections/proj/tpeqd.hpp

Tags: upstream-0.15.3+svn20934
Import upstream version 0.15.3+svn20934

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/builtin-ggl/ggl/extensions/gis/projections/proj/tpeqd.hpp
 
1
#ifndef GGL_PROJECTIONS_TPEQD_HPP
 
2
#define GGL_PROJECTIONS_TPEQD_HPP
 
3
 
 
4
// Generic Geometry Library - projections (based on PROJ4)
 
5
// This file is automatically generated. DO NOT EDIT.
 
6
 
 
7
// Copyright Barend Gehrels (1995-2009), Geodan Holding B.V. Amsterdam, the Netherlands.
 
8
// Copyright Bruno Lalande (2008-2009)
 
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 Geometry Library by Barend Gehrels (Geodan, Amsterdam)
 
17
 
 
18
// Original copyright notice:
 
19
 
 
20
// Permission is hereby granted, free of charge, to any person obtaining a
 
21
// copy of this software and associated documentation files (the "Software"),
 
22
// to deal in the Software without restriction, including without limitation
 
23
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
 
24
// and/or sell copies of the Software, and to permit persons to whom the
 
25
// Software is furnished to do so, subject to the following conditions:
 
26
 
 
27
// The above copyright notice and this permission notice shall be included
 
28
// in all copies or substantial portions of the Software.
 
29
 
 
30
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 
31
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
32
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 
33
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
34
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
35
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 
36
// DEALINGS IN THE SOFTWARE.
 
37
 
 
38
#include <boost/math/special_functions/hypot.hpp>
 
39
 
 
40
#include <ggl/extensions/gis/projections/impl/base_static.hpp>
 
41
#include <ggl/extensions/gis/projections/impl/base_dynamic.hpp>
 
42
#include <ggl/extensions/gis/projections/impl/projects.hpp>
 
43
#include <ggl/extensions/gis/projections/impl/factory_entry.hpp>
 
44
 
 
45
namespace ggl { namespace projection
 
46
{
 
47
    #ifndef DOXYGEN_NO_DETAIL
 
48
    namespace detail { namespace tpeqd{ 
 
49
 
 
50
            struct par_tpeqd
 
51
            {
 
52
                double cp1, sp1, cp2, sp2, ccs, cs, sc, r2z0, z02, dlam2;
 
53
                double hz0, thz0, rhshz0, ca, sa, lp, lamc;
 
54
            };
 
55
 
 
56
            // template class, using CRTP to implement forward/inverse
 
57
            template <typename Geographic, typename Cartesian, typename Parameters>
 
58
            struct base_tpeqd_spheroid : public base_t_fi<base_tpeqd_spheroid<Geographic, Cartesian, Parameters>,
 
59
                     Geographic, Cartesian, Parameters>
 
60
            {
 
61
 
 
62
                 typedef double geographic_type;
 
63
                 typedef double cartesian_type;
 
64
 
 
65
                par_tpeqd m_proj_parm;
 
66
 
 
67
                inline base_tpeqd_spheroid(const Parameters& par)
 
68
                    : base_t_fi<base_tpeqd_spheroid<Geographic, Cartesian, Parameters>,
 
69
                     Geographic, Cartesian, Parameters>(*this, par) {}
 
70
 
 
71
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 
72
                {
 
73
                    double t, z1, z2, dl1, dl2, sp, cp;
 
74
                
 
75
                    sp = sin(lp_lat);
 
76
                    cp = cos(lp_lat);
 
77
                    z1 = aacos(this->m_proj_parm.sp1 * sp + this->m_proj_parm.cp1 * cp * cos(dl1 = lp_lon + this->m_proj_parm.dlam2));
 
78
                    z2 = aacos(this->m_proj_parm.sp2 * sp + this->m_proj_parm.cp2 * cp * cos(dl2 = lp_lon - this->m_proj_parm.dlam2));
 
79
                    z1 *= z1;
 
80
                    z2 *= z2;
 
81
                    xy_x = this->m_proj_parm.r2z0 * (t = z1 - z2);
 
82
                    t = this->m_proj_parm.z02 - t;
 
83
                    xy_y = this->m_proj_parm.r2z0 * asqrt(4. * this->m_proj_parm.z02 * z2 - t * t);
 
84
                    if ((this->m_proj_parm.ccs * sp - cp * (this->m_proj_parm.cs * sin(dl1) - this->m_proj_parm.sc * sin(dl2))) < 0.)
 
85
                        xy_y = -xy_y;
 
86
                }
 
87
 
 
88
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 
89
                {
 
90
                    double cz1, cz2, s, d, cp, sp;
 
91
                
 
92
                    cz1 = cos(boost::math::hypot(xy_y, xy_x + this->m_proj_parm.hz0));
 
93
                    cz2 = cos(boost::math::hypot(xy_y, xy_x - this->m_proj_parm.hz0));
 
94
                    s = cz1 + cz2;
 
95
                    d = cz1 - cz2;
 
96
                    lp_lon = - atan2(d, (s * this->m_proj_parm.thz0));
 
97
                    lp_lat = aacos(boost::math::hypot(this->m_proj_parm.thz0 * s, d) * this->m_proj_parm.rhshz0);
 
98
                    if ( xy_y < 0. )
 
99
                        lp_lat = - lp_lat;
 
100
                    /* lam--phi now in system relative to P1--P2 base equator */
 
101
                    sp = sin(lp_lat);
 
102
                    cp = cos(lp_lat);
 
103
                    lp_lat = aasin(this->m_proj_parm.sa * sp + this->m_proj_parm.ca * cp * (s = cos(lp_lon -= this->m_proj_parm.lp)));
 
104
                    lp_lon = atan2(cp * sin(lp_lon), this->m_proj_parm.sa * cp * s - this->m_proj_parm.ca * sp) + this->m_proj_parm.lamc;
 
105
                }
 
106
            };
 
107
 
 
108
            // Two Point Equidistant
 
109
            template <typename Parameters>
 
110
            void setup_tpeqd(Parameters& par, par_tpeqd& proj_parm)
 
111
            {
 
112
                double lam_1, lam_2, phi_1, phi_2, A12, pp;
 
113
                /* get control point locations */
 
114
                phi_1 = pj_param(par.params, "rlat_1").f;
 
115
                lam_1 = pj_param(par.params, "rlon_1").f;
 
116
                phi_2 = pj_param(par.params, "rlat_2").f;
 
117
                lam_2 = pj_param(par.params, "rlon_2").f;
 
118
                if (phi_1 == phi_2 && lam_1 == lam_2) throw proj_exception(-25);
 
119
                par.lam0 = adjlon(0.5 * (lam_1 + lam_2));
 
120
                proj_parm.dlam2 = adjlon(lam_2 - lam_1);
 
121
                proj_parm.cp1 = cos(phi_1);
 
122
                proj_parm.cp2 = cos(phi_2);
 
123
                proj_parm.sp1 = sin(phi_1);
 
124
                proj_parm.sp2 = sin(phi_2);
 
125
                proj_parm.cs = proj_parm.cp1 * proj_parm.sp2;
 
126
                proj_parm.sc = proj_parm.sp1 * proj_parm.cp2;
 
127
                proj_parm.ccs = proj_parm.cp1 * proj_parm.cp2 * sin(proj_parm.dlam2);
 
128
                proj_parm.z02 = aacos(proj_parm.sp1 * proj_parm.sp2 + proj_parm.cp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
 
129
                proj_parm.hz0 = .5 * proj_parm.z02;
 
130
                A12 = atan2(proj_parm.cp2 * sin(proj_parm.dlam2),
 
131
                    proj_parm.cp1 * proj_parm.sp2 - proj_parm.sp1 * proj_parm.cp2 * cos(proj_parm.dlam2));
 
132
                proj_parm.ca = cos(pp = aasin(proj_parm.cp1 * sin(A12)));
 
133
                proj_parm.sa = sin(pp);
 
134
                proj_parm.lp = adjlon(atan2(proj_parm.cp1 * cos(A12), proj_parm.sp1) - proj_parm.hz0);
 
135
                proj_parm.dlam2 *= .5;
 
136
                proj_parm.lamc = HALFPI - atan2(sin(A12) * proj_parm.sp1, cos(A12)) - proj_parm.dlam2;
 
137
                proj_parm.thz0 = tan(proj_parm.hz0);
 
138
                proj_parm.rhshz0 = .5 / sin(proj_parm.hz0);
 
139
                proj_parm.r2z0 = 0.5 / proj_parm.z02;
 
140
                proj_parm.z02 *= proj_parm.z02;
 
141
                // par.inv = s_inverse;
 
142
                // par.fwd = s_forward;
 
143
                par.es = 0.;
 
144
            }
 
145
 
 
146
        }} // namespace detail::tpeqd
 
147
    #endif // doxygen 
 
148
 
 
149
    /*!
 
150
        \brief Two Point Equidistant projection
 
151
        \ingroup projections
 
152
        \tparam Geographic latlong point type
 
153
        \tparam Cartesian xy point type
 
154
        \tparam Parameters parameter type
 
155
        \par Projection characteristics
 
156
         - Miscellaneous
 
157
         - Spheroid
 
158
         - lat_1= lon_1= lat_2= lon_2=
 
159
        \par Example
 
160
        \image html ex_tpeqd.gif
 
161
    */
 
162
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
 
163
    struct tpeqd_spheroid : public detail::tpeqd::base_tpeqd_spheroid<Geographic, Cartesian, Parameters>
 
164
    {
 
165
        inline tpeqd_spheroid(const Parameters& par) : detail::tpeqd::base_tpeqd_spheroid<Geographic, Cartesian, Parameters>(par)
 
166
        {
 
167
            detail::tpeqd::setup_tpeqd(this->m_par, this->m_proj_parm);
 
168
        }
 
169
    };
 
170
 
 
171
    #ifndef DOXYGEN_NO_DETAIL
 
172
    namespace detail
 
173
    {
 
174
 
 
175
        // Factory entry(s)
 
176
        template <typename Geographic, typename Cartesian, typename Parameters>
 
177
        class tpeqd_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
 
178
        {
 
179
            public :
 
180
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
 
181
                {
 
182
                    return new base_v_fi<tpeqd_spheroid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
 
183
                }
 
184
        };
 
185
 
 
186
        template <typename Geographic, typename Cartesian, typename Parameters>
 
187
        inline void tpeqd_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
 
188
        {
 
189
            factory.add_to_factory("tpeqd", new tpeqd_entry<Geographic, Cartesian, Parameters>);
 
190
        }
 
191
 
 
192
    } // namespace detail 
 
193
    #endif // doxygen
 
194
 
 
195
}} // namespace ggl::projection
 
196
 
 
197
#endif // GGL_PROJECTIONS_TPEQD_HPP
 
198