← Back to branch summary

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

« back to all changes in this revision

Viewing changes to include/ggl/projections/proj/imw_p.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Bernd Zeimetz
  • Date: 2009-09-13 00:52:12 UTC
  • mto: (1.2.7 upstream) (0.1.3 upstream) (3.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090913005212-pjecal8zxm07x0fj
Import upstream version 0.14+svnfixes~20090912

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/ggl/projections/proj/imw_p.hpp
 
1
#ifndef GGL_PROJECTIONS_IMW_P_HPP
 
2
#define GGL_PROJECTIONS_IMW_P_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/projections/impl/base_static.hpp>
 
41
#include <ggl/projections/impl/base_dynamic.hpp>
 
42
#include <ggl/projections/impl/projects.hpp>
 
43
#include <ggl/projections/impl/factory_entry.hpp>
 
44
#include <ggl/projections/impl/pj_mlfn.hpp>
 
45
 
 
46
namespace ggl { namespace projection
 
47
{
 
48
    #ifndef DOXYGEN_NO_DETAIL
 
49
    namespace detail { namespace imw_p{ 
 
50
            static const double TOL = 1e-10;
 
51
            static const double EPS = 1e-10;
 
52
 
 
53
            struct PXY { double x, y; }; // x/y projection specific
 
54
 
 
55
            struct par_imw_p
 
56
            {
 
57
                double  P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
 
58
                double  phi_1, phi_2, lam_1;
 
59
                double en[EN_SIZE];
 
60
                int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
 
61
            };
 
62
            template <typename Parameters>
 
63
                inline int
 
64
            phi12(Parameters& par, par_imw_p& proj_parm, double *del, double *sig) {
 
65
                int err = 0;
 
66
            
 
67
                if (!pj_param(par.params, "tlat_1").i ||
 
68
                    !pj_param(par.params, "tlat_2").i) {
 
69
                    err = -41;
 
70
                } else {
 
71
                    proj_parm.phi_1 = pj_param(par.params, "rlat_1").f;
 
72
                    proj_parm.phi_2 = pj_param(par.params, "rlat_2").f;
 
73
                    *del = 0.5 * (proj_parm.phi_2 - proj_parm.phi_1);
 
74
                    *sig = 0.5 * (proj_parm.phi_2 + proj_parm.phi_1);
 
75
                    err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
 
76
                }
 
77
                return err;
 
78
            }
 
79
            template <typename Parameters>
 
80
                inline PXY
 
81
            loc_for(const double& lp_lam, const double& lp_phi, const Parameters& par, const par_imw_p& proj_parm, double *yc) {
 
82
                PXY xy;
 
83
            
 
84
                if (! lp_phi) {
 
85
                    xy.x = lp_lam;
 
86
                    xy.y = 0.;
 
87
                } else {
 
88
                    double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
 
89
            
 
90
                    sp = sin(lp_phi);
 
91
                    m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
 
92
                    xa = proj_parm.Pp + proj_parm.Qp * m;
 
93
                    ya = proj_parm.P + proj_parm.Q * m;
 
94
                    R = 1. / (tan(lp_phi) * sqrt(1. - par.es * sp * sp));
 
95
                    C = sqrt(R * R - xa * xa);
 
96
                    if (lp_phi < 0.) C = - C;
 
97
                    C += ya - R;
 
98
                    if (proj_parm.mode < 0) {
 
99
                        xb = lp_lam;
 
100
                        yb = proj_parm.C2;
 
101
                    } else {
 
102
                        t = lp_lam * proj_parm.sphi_2;
 
103
                        xb = proj_parm.R_2 * sin(t);
 
104
                        yb = proj_parm.C2 + proj_parm.R_2 * (1. - cos(t));
 
105
                    }
 
106
                    if (proj_parm.mode > 0) {
 
107
                        xc = lp_lam;
 
108
                        *yc = 0.;
 
109
                    } else {
 
110
                        t = lp_lam * proj_parm.sphi_1;
 
111
                        xc = proj_parm.R_1 * sin(t);
 
112
                        *yc = proj_parm.R_1 * (1. - cos(t));
 
113
                    }
 
114
                    D = (xb - xc)/(yb - *yc);
 
115
                    B = xc + D * (C + R - *yc);
 
116
                    xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
 
117
                    if (lp_phi > 0)
 
118
                        xy.x = - xy.x;
 
119
                    xy.x = (B + xy.x) / (1. + D * D);
 
120
                    xy.y = sqrt(R * R - xy.x * xy.x);
 
121
                    if (lp_phi > 0)
 
122
                        xy.y = - xy.y;
 
123
                    xy.y += C + R;
 
124
                }
 
125
                return (xy);
 
126
            }
 
127
            
 
128
            template <typename Parameters>
 
129
            inline void
 
130
            xy(Parameters& par, par_imw_p& proj_parm, double phi, double *x, double *y, double *sp, double *R) {
 
131
                double F;
 
132
            
 
133
                *sp = sin(phi);
 
134
                *R = 1./(tan(phi) * sqrt(1. - par.es * *sp * *sp ));
 
135
                F = proj_parm.lam_1 * *sp;
 
136
                *y = *R * (1 - cos(F));
 
137
                *x = *R * sin(F);
 
138
            }
 
139
            
 
140
 
 
141
            // template class, using CRTP to implement forward/inverse
 
142
            template <typename Geographic, typename Cartesian, typename Parameters>
 
143
            struct base_imw_p_ellipsoid : public base_t_fi<base_imw_p_ellipsoid<Geographic, Cartesian, Parameters>,
 
144
                     Geographic, Cartesian, Parameters>
 
145
            {
 
146
 
 
147
                 typedef double geographic_type;
 
148
                 typedef double cartesian_type;
 
149
 
 
150
                par_imw_p m_proj_parm;
 
151
 
 
152
                inline base_imw_p_ellipsoid(const Parameters& par)
 
153
                    : base_t_fi<base_imw_p_ellipsoid<Geographic, Cartesian, Parameters>,
 
154
                     Geographic, Cartesian, Parameters>(*this, par) {}
 
155
 
 
156
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 
157
                {
 
158
                    double yc = 0;
 
159
                    PXY xy = loc_for(lp_lon, lp_lat, this->m_par, m_proj_parm, &yc);
 
160
                        xy_x = xy.x; xy_y = xy.y;
 
161
                }
 
162
 
 
163
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 
164
                {
 
165
                    PXY t;
 
166
                    double yc = 0;
 
167
                
 
168
                    lp_lat = this->m_proj_parm.phi_2;
 
169
                    lp_lon = xy_x / cos(lp_lat);
 
170
                    do {
 
171
                        t = loc_for(lp_lon, lp_lat, this->m_par, m_proj_parm, &yc);
 
172
                        lp_lat = ((lp_lat - this->m_proj_parm.phi_1) * (xy_y - yc) / (t.y - yc)) + this->m_proj_parm.phi_1;
 
173
                        lp_lon = lp_lon * xy_x / t.x;
 
174
                    } while (fabs(t.x - xy_x) > TOL || fabs(t.y - xy_y) > TOL);
 
175
                }
 
176
            };
 
177
 
 
178
            // International Map of the World Polyconic
 
179
            template <typename Parameters>
 
180
            void setup_imw_p(Parameters& par, par_imw_p& proj_parm)
 
181
            {
 
182
                double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
 
183
                int i;
 
184
                    pj_enfn(par.es, proj_parm.en);
 
185
                if( (i = phi12(par, proj_parm, &del, &sig)) != 0)
 
186
                    throw proj_exception(i);
 
187
                if (proj_parm.phi_2 < proj_parm.phi_1) { /* make sure proj_parm.phi_1 most southerly */
 
188
                    del = proj_parm.phi_1;
 
189
                    proj_parm.phi_1 = proj_parm.phi_2;
 
190
                    proj_parm.phi_2 = del;
 
191
                }
 
192
                if (pj_param(par.params, "tlon_1").i)
 
193
                    proj_parm.lam_1 = pj_param(par.params, "rlon_1").f;
 
194
                else { /* use predefined based upon latitude */
 
195
                    sig = fabs(sig * RAD_TO_DEG);
 
196
                    if (sig <= 60)      sig = 2.;
 
197
                    else if (sig <= 76) sig = 4.;
 
198
                    else                sig = 8.;
 
199
                    proj_parm.lam_1 = sig * DEG_TO_RAD;
 
200
                }
 
201
                proj_parm.mode = 0;
 
202
                if (proj_parm.phi_1) xy(par, proj_parm, proj_parm.phi_1, &x1, &y1, &proj_parm.sphi_1, &proj_parm.R_1);
 
203
                else {
 
204
                    proj_parm.mode = 1;
 
205
                    y1 = 0.;
 
206
                    x1 = proj_parm.lam_1;
 
207
                }
 
208
                if (proj_parm.phi_2) xy(par, proj_parm, proj_parm.phi_2, &x2, &T2, &proj_parm.sphi_2, &proj_parm.R_2);
 
209
                else {
 
210
                    proj_parm.mode = -1;
 
211
                    T2 = 0.;
 
212
                    x2 = proj_parm.lam_1;
 
213
                }
 
214
                m1 = pj_mlfn(proj_parm.phi_1, proj_parm.sphi_1, cos(proj_parm.phi_1), proj_parm.en);
 
215
                m2 = pj_mlfn(proj_parm.phi_2, proj_parm.sphi_2, cos(proj_parm.phi_2), proj_parm.en);
 
216
                t = m2 - m1;
 
217
                s = x2 - x1;
 
218
                y2 = sqrt(t * t - s * s) + y1;
 
219
                proj_parm.C2 = y2 - T2;
 
220
                t = 1. / t;
 
221
                proj_parm.P = (m2 * y1 - m1 * y2) * t;
 
222
                proj_parm.Q = (y2 - y1) * t;
 
223
                proj_parm.Pp = (m2 * x1 - m1 * x2) * t;
 
224
                proj_parm.Qp = (x2 - x1) * t;
 
225
                // par.fwd = e_forward;
 
226
                // par.inv = e_inverse;
 
227
            }
 
228
 
 
229
        }} // namespace detail::imw_p
 
230
    #endif // doxygen 
 
231
 
 
232
    /*!
 
233
        \brief International Map of the World Polyconic projection
 
234
        \ingroup projections
 
235
        \tparam Geographic latlong point type
 
236
        \tparam Cartesian xy point type
 
237
        \tparam Parameters parameter type
 
238
        \par Projection characteristics
 
239
         - Mod Polyconic
 
240
         - Ellipsoid
 
241
         - lat_1= and lat_2= [lon_1=]
 
242
        \par Example
 
243
        \image html ex_imw_p.gif
 
244
    */
 
245
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
 
246
    struct imw_p_ellipsoid : public detail::imw_p::base_imw_p_ellipsoid<Geographic, Cartesian, Parameters>
 
247
    {
 
248
        inline imw_p_ellipsoid(const Parameters& par) : detail::imw_p::base_imw_p_ellipsoid<Geographic, Cartesian, Parameters>(par)
 
249
        {
 
250
            detail::imw_p::setup_imw_p(this->m_par, this->m_proj_parm);
 
251
        }
 
252
    };
 
253
 
 
254
    #ifndef DOXYGEN_NO_DETAIL
 
255
    namespace detail
 
256
    {
 
257
 
 
258
        // Factory entry(s)
 
259
        template <typename Geographic, typename Cartesian, typename Parameters>
 
260
        class imw_p_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
 
261
        {
 
262
            public :
 
263
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
 
264
                {
 
265
                    return new base_v_fi<imw_p_ellipsoid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
 
266
                }
 
267
        };
 
268
 
 
269
        template <typename Geographic, typename Cartesian, typename Parameters>
 
270
        inline void imw_p_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
 
271
        {
 
272
            factory.add_to_factory("imw_p", new imw_p_entry<Geographic, Cartesian, Parameters>);
 
273
        }
 
274
 
 
275
    } // namespace detail 
 
276
    #endif // doxygen
 
277
 
 
278
}} // namespace ggl::projection
 
279
 
 
280
#endif // GGL_PROJECTIONS_IMW_P_HPP
 
281