← Back to branch summary

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

« back to all changes in this revision

Viewing changes to include/ggl/projections/proj/bipc.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/bipc.hpp
 
1
#ifndef GGL_PROJECTIONS_BIPC_HPP
 
2
#define GGL_PROJECTIONS_BIPC_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
 
 
45
namespace ggl { namespace projection
 
46
{
 
47
    #ifndef DOXYGEN_NO_DETAIL
 
48
    namespace detail { namespace bipc{ 
 
49
            static const double EPS = 1e-10;
 
50
            static const double EPS10 = 1e-10;
 
51
            static const double ONEEPS = 1.000000001;
 
52
            static const int NITER = 10;
 
53
            static const double lamB = -.34894976726250681539;
 
54
            static const double n = .63055844881274687180;
 
55
            static const double F = 1.89724742567461030582;
 
56
            static const double Azab = .81650043674686363166;
 
57
            static const double Azba = 1.82261843856185925133;
 
58
            static const double T = 1.27246578267089012270;
 
59
            static const double rhoc = 1.20709121521568721927;
 
60
            static const double cAzc = .69691523038678375519;
 
61
            static const double sAzc = .71715351331143607555;
 
62
            static const double C45 = .70710678118654752469;
 
63
            static const double S45 = .70710678118654752410;
 
64
            static const double C20 = .93969262078590838411;
 
65
            static const double S20 = -.34202014332566873287;
 
66
            static const double R110 = 1.91986217719376253360;
 
67
            static const double R104 = 1.81514242207410275904;
 
68
 
 
69
            struct par_bipc
 
70
            {
 
71
                int noskew;
 
72
            };
 
73
 
 
74
            // template class, using CRTP to implement forward/inverse
 
75
            template <typename Geographic, typename Cartesian, typename Parameters>
 
76
            struct base_bipc_spheroid : public base_t_fi<base_bipc_spheroid<Geographic, Cartesian, Parameters>,
 
77
                     Geographic, Cartesian, Parameters>
 
78
            {
 
79
 
 
80
                 typedef double geographic_type;
 
81
                 typedef double cartesian_type;
 
82
 
 
83
                par_bipc m_proj_parm;
 
84
 
 
85
                inline base_bipc_spheroid(const Parameters& par)
 
86
                    : base_t_fi<base_bipc_spheroid<Geographic, Cartesian, Parameters>,
 
87
                     Geographic, Cartesian, Parameters>(*this, par) {}
 
88
 
 
89
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 
90
                {
 
91
                        double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
 
92
                        int tag;
 
93
                
 
94
                        cphi = cos(lp_lat);
 
95
                        sphi = sin(lp_lat);
 
96
                        cdlam = cos(sdlam = lamB - lp_lon);
 
97
                        sdlam = sin(sdlam);
 
98
                        if (fabs(fabs(lp_lat) - HALFPI) < EPS10) {
 
99
                                Az = lp_lat < 0. ? PI : 0.;
 
100
                                tphi = HUGE_VAL;
 
101
                        } else {
 
102
                                tphi = sphi / cphi;
 
103
                                Az = atan2(sdlam , C45 * (tphi - cdlam));
 
104
                        }
 
105
                        if( (tag = (Az > Azba)) ) {
 
106
                                cdlam = cos(sdlam = lp_lon + R110);
 
107
                                sdlam = sin(sdlam);
 
108
                                z = S20 * sphi + C20 * cphi * cdlam;
 
109
                                if (fabs(z) > 1.) {
 
110
                                        if (fabs(z) > ONEEPS) throw proj_exception();
 
111
                                        else z = z < 0. ? -1. : 1.;
 
112
                                } else
 
113
                                        z = acos(z);
 
114
                                if (tphi != HUGE_VAL)
 
115
                                        Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
 
116
                                Av = Azab;
 
117
                                xy_y = rhoc;
 
118
                        } else {
 
119
                                z = S45 * (sphi + cphi * cdlam);
 
120
                                if (fabs(z) > 1.) {
 
121
                                        if (fabs(z) > ONEEPS) throw proj_exception();
 
122
                                        else z = z < 0. ? -1. : 1.;
 
123
                                } else
 
124
                                        z = acos(z);
 
125
                                Av = Azba;
 
126
                                xy_y = -rhoc;
 
127
                        }
 
128
                        if (z < 0.) throw proj_exception();;
 
129
                        r = F * (t = pow(tan(.5 * z), n));
 
130
                        if ((al = .5 * (R104 - z)) < 0.) throw proj_exception();;
 
131
                        al = (t + pow(al, n)) / T;
 
132
                        if (fabs(al) > 1.) {
 
133
                                if (fabs(al) > ONEEPS) throw proj_exception();
 
134
                                else al = al < 0. ? -1. : 1.;
 
135
                        } else
 
136
                                al = acos(al);
 
137
                        if (fabs(t = n * (Av - Az)) < al)
 
138
                                r /= cos(al + (tag ? t : -t));
 
139
                        xy_x = r * sin(t);
 
140
                        xy_y += (tag ? -r : r) * cos(t);
 
141
                        if (this->m_proj_parm.noskew) {
 
142
                                t = xy_x;
 
143
                                xy_x = -xy_x * cAzc - xy_y * sAzc; 
 
144
                                xy_y = -xy_y * cAzc + t * sAzc; 
 
145
                        }
 
146
                }
 
147
 
 
148
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 
149
                {
 
150
                        double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
 
151
                        int neg, i;
 
152
                
 
153
                        if (this->m_proj_parm.noskew) {
 
154
                                t = xy_x;
 
155
                                xy_x = -xy_x * cAzc + xy_y * sAzc; 
 
156
                                xy_y = -xy_y * cAzc - t * sAzc; 
 
157
                        }
 
158
                        if( (neg = (xy_x < 0.)) ) {
 
159
                                xy_y = rhoc - xy_y;
 
160
                                s = S20;
 
161
                                c = C20;
 
162
                                Av = Azab;
 
163
                        } else {
 
164
                                xy_y += rhoc;
 
165
                                s = S45;
 
166
                                c = C45;
 
167
                                Av = Azba;
 
168
                        }
 
169
                        rl = rp = r = boost::math::hypot(xy_x, xy_y);
 
170
                        fAz = fabs(Az = atan2(xy_x, xy_y));
 
171
                        for (i = NITER; i ; --i) {
 
172
                                z = 2. * atan(pow(r / F,1 / n));
 
173
                                al = acos((pow(tan(.5 * z), n) +
 
174
                                   pow(tan(.5 * (R104 - z)), n)) / T);
 
175
                                if (fAz < al)
 
176
                                        r = rp * cos(al + (neg ? Az : -Az));
 
177
                                if (fabs(rl - r) < EPS)
 
178
                                        break;
 
179
                                rl = r;
 
180
                        }
 
181
                        if (! i) throw proj_exception();;
 
182
                        Az = Av - Az / n;
 
183
                        lp_lat = asin(s * cos(z) + c * sin(z) * cos(Az));
 
184
                        lp_lon = atan2(sin(Az), c / tan(z) - s * cos(Az));
 
185
                        if (neg)
 
186
                                lp_lon -= R110;
 
187
                        else
 
188
                                lp_lon = lamB - lp_lon;
 
189
                }
 
190
            };
 
191
 
 
192
            // Bipolar conic of western hemisphere
 
193
            template <typename Parameters>
 
194
            void setup_bipc(Parameters& par, par_bipc& proj_parm)
 
195
            {
 
196
                proj_parm.noskew = pj_param(par.params, "bns").i;
 
197
                // par.inv = s_inverse;
 
198
                // par.fwd = s_forward;
 
199
                par.es = 0.;
 
200
            }
 
201
 
 
202
        }} // namespace detail::bipc
 
203
    #endif // doxygen 
 
204
 
 
205
    /*!
 
206
        \brief Bipolar conic of western hemisphere projection
 
207
        \ingroup projections
 
208
        \tparam Geographic latlong point type
 
209
        \tparam Cartesian xy point type
 
210
        \tparam Parameters parameter type
 
211
        \par Projection characteristics
 
212
         - Conic
 
213
         - Spheroid
 
214
        \par Example
 
215
        \image html ex_bipc.gif
 
216
    */
 
217
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
 
218
    struct bipc_spheroid : public detail::bipc::base_bipc_spheroid<Geographic, Cartesian, Parameters>
 
219
    {
 
220
        inline bipc_spheroid(const Parameters& par) : detail::bipc::base_bipc_spheroid<Geographic, Cartesian, Parameters>(par)
 
221
        {
 
222
            detail::bipc::setup_bipc(this->m_par, this->m_proj_parm);
 
223
        }
 
224
    };
 
225
 
 
226
    #ifndef DOXYGEN_NO_DETAIL
 
227
    namespace detail
 
228
    {
 
229
 
 
230
        // Factory entry(s)
 
231
        template <typename Geographic, typename Cartesian, typename Parameters>
 
232
        class bipc_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
 
233
        {
 
234
            public :
 
235
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
 
236
                {
 
237
                    return new base_v_fi<bipc_spheroid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
 
238
                }
 
239
        };
 
240
 
 
241
        template <typename Geographic, typename Cartesian, typename Parameters>
 
242
        inline void bipc_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
 
243
        {
 
244
            factory.add_to_factory("bipc", new bipc_entry<Geographic, Cartesian, Parameters>);
 
245
        }
 
246
 
 
247
    } // namespace detail 
 
248
    #endif // doxygen
 
249
 
 
250
}} // namespace ggl::projection
 
251
 
 
252
#endif // GGL_PROJECTIONS_BIPC_HPP
 
253