Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:35:43 

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002 
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004 
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008 
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012 
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017 
0018 // Last updated version of proj: 5.0.0
0019 
0020 // Original copyright notice:
0021 
0022 // This code was entirely written by Nathan Wagner
0023 // and is in the public domain.
0024 
0025 // Permission is hereby granted, free of charge, to any person obtaining a
0026 // copy of this software and associated documentation files (the "Software"),
0027 // to deal in the Software without restriction, including without limitation
0028 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0029 // and/or sell copies of the Software, and to permit persons to whom the
0030 // Software is furnished to do so, subject to the following conditions:
0031 
0032 // The above copyright notice and this permission notice shall be included
0033 // in all copies or substantial portions of the Software.
0034 
0035 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0036 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0037 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0038 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0039 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0040 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0041 // DEALINGS IN THE SOFTWARE.
0042 
0043 #ifndef BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
0044 #define BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
0045 
0046 #include <sstream>
0047 
0048 #include <boost/core/ignore_unused.hpp>
0049 
0050 #include <boost/geometry/core/assert.hpp>
0051 
0052 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0053 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0054 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0055 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0056 #include <boost/geometry/srs/projections/impl/projects.hpp>
0057 
0058 #include <boost/geometry/util/math.hpp>
0059 
0060 namespace boost { namespace geometry
0061 {
0062 
0063 namespace projections
0064 {
0065     #ifndef DOXYGEN_NO_DETAIL
0066     namespace detail { namespace isea
0067     {
0068             static const double epsilon = std::numeric_limits<double>::epsilon();
0069 
0070             /* sqrt(5)/M_PI */
0071             static const double isea_scale = 0.8301572857837594396028083;
0072             /* 26.565051177 degrees */
0073             static const double v_lat = 0.46364760899944494524;
0074             /* 52.62263186 */
0075             static const double e_rad = 0.91843818702186776133;
0076             /* 10.81231696 */
0077             static const double f_rad = 0.18871053072122403508;
0078             /* R tan(g) sin(60) */
0079             static const double table_g = 0.6615845383;
0080             /* H = 0.25 R tan g = */
0081             static const double table_h = 0.1909830056;
0082             //static const double RPRIME = 0.91038328153090290025;
0083             static const double isea_std_lat = 1.01722196792335072101;
0084             static const double isea_std_lon = .19634954084936207740;
0085 
0086             template <typename T>
0087             inline T deg30_rad() { return T(30) * geometry::math::d2r<T>(); }
0088             template <typename T>
0089             inline T deg120_rad() { return T(120) * geometry::math::d2r<T>(); }
0090             template <typename T>
0091             inline T deg72_rad() { return T(72) * geometry::math::d2r<T>(); }
0092             template <typename T>
0093             inline T deg90_rad() { return geometry::math::half_pi<T>(); }
0094             template <typename T>
0095             inline T deg144_rad() { return T(144) * geometry::math::d2r<T>(); }
0096             template <typename T>
0097             inline T deg36_rad() { return T(36) * geometry::math::d2r<T>(); }
0098             template <typename T>
0099             inline T deg108_rad() { return T(108) * geometry::math::d2r<T>(); }
0100             template <typename T>
0101             inline T deg180_rad() { return geometry::math::pi<T>(); }
0102 
0103             inline bool downtri(int tri) { return (((tri - 1) / 5) % 2 == 1); }
0104 
0105             /*
0106              * Proj 4 provides its own entry points into
0107              * the code, so none of the library functions
0108              * need to be global
0109              */
0110 
0111             struct hex {
0112                     int iso;
0113                     int x, y, z;
0114             };
0115 
0116             /* y *must* be positive down as the xy /iso conversion assumes this */
0117             inline
0118             int hex_xy(struct hex *h) {
0119                 if (!h->iso) return 1;
0120                 if (h->x >= 0) {
0121                     h->y = -h->y - (h->x+1)/2;
0122                 } else {
0123                     /* need to round toward -inf, not toward zero, so x-1 */
0124                     h->y = -h->y - h->x/2;
0125                 }
0126                 h->iso = 0;
0127 
0128                 return 1;
0129             }
0130 
0131             inline
0132             int hex_iso(struct hex *h) {
0133                 if (h->iso) return 1;
0134 
0135                 if (h->x >= 0) {
0136                     h->y = (-h->y - (h->x+1)/2);
0137                 } else {
0138                     /* need to round toward -inf, not toward zero, so x-1 */
0139                     h->y = (-h->y - (h->x)/2);
0140                 }
0141 
0142                 h->z = -h->x - h->y;
0143                 h->iso = 1;
0144                 return 1;
0145             }
0146 
0147             template <typename T>
0148             inline
0149             int hexbin2(T const& width, T x, T y, int *i, int *j)
0150             {
0151                 T z, rx, ry, rz;
0152                 T abs_dx, abs_dy, abs_dz;
0153                 int ix, iy, iz, s;
0154                 struct hex h;
0155 
0156                 static const T cos_deg30 = cos(deg30_rad<T>());
0157 
0158                 x = x / cos_deg30; /* rotated X coord */
0159                 y = y - x / 2.0; /* adjustment for rotated X */
0160 
0161                 /* adjust for actual hexwidth */
0162                 x /= width;
0163                 y /= width;
0164 
0165                 z = -x - y;
0166 
0167                 rx = floor(x + 0.5);
0168                 ix = (int)rx;
0169                 ry = floor(y + 0.5);
0170                 iy = (int)ry;
0171                 rz = floor(z + 0.5);
0172                 iz = (int)rz;
0173 
0174                 s = ix + iy + iz;
0175 
0176                 if (s) {
0177                     abs_dx = fabs(rx - x);
0178                     abs_dy = fabs(ry - y);
0179                     abs_dz = fabs(rz - z);
0180 
0181                     if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
0182                         ix -= s;
0183                     } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
0184                         iy -= s;
0185                     } else {
0186                         iz -= s;
0187                     }
0188                 }
0189                 h.x = ix;
0190                 h.y = iy;
0191                 h.z = iz;
0192                 h.iso = 1;
0193 
0194                 hex_xy(&h);
0195                 *i = h.x;
0196                 *j = h.y;
0197                 return ix * 100 + iy;
0198             }
0199 
0200             //enum isea_poly { isea_none = 0, isea_icosahedron = 20 };
0201             //enum isea_topology { isea_hexagon=6, isea_triangle=3, isea_diamond=4 };
0202             enum isea_address_form {
0203                 /*isea_addr_geo,*/ isea_addr_q2di, isea_addr_seqnum,
0204                 /*isea_addr_interleave,*/ isea_addr_plane, isea_addr_q2dd,
0205                 isea_addr_projtri, isea_addr_vertex2dd, isea_addr_hex
0206             };
0207 
0208             template <typename T>
0209             struct isea_dgg {
0210                 T                 o_lat, o_lon, o_az; /* orientation, radians */
0211                 T                 radius; /* radius of the earth in meters, ignored 1.0 */
0212                 unsigned long     serial;
0213                 //int               pole; /* true if standard snyder */
0214                 int               aperture; /* valid values depend on partitioning method */
0215                 int               resolution;
0216                 int               triangle; /* triangle of last transformed point */
0217                 int               quad; /* quad of last transformed point */
0218                 //isea_poly         polyhedron; /* ignored, icosahedron */
0219                 //isea_topology     topology; /* ignored, hexagon */
0220                 isea_address_form output; /* an isea_address_form */
0221             };
0222 
0223             template <typename T>
0224             struct isea_pt {
0225                 T x, y;
0226             };
0227 
0228             template <typename T>
0229             struct isea_geo {
0230                 T lon, lat;
0231             };
0232 
0233             template <typename T>
0234             struct isea_address {
0235                 T      x,y; /* or i,j or lon,lat depending on type */
0236                 int    type; /* enum isea_address_form */
0237                 int    number;
0238             };
0239 
0240             /* ENDINC */
0241 
0242             enum snyder_polyhedron {
0243                 snyder_poly_hexagon = 0, snyder_poly_pentagon = 1,
0244                 snyder_poly_tetrahedron = 2, snyder_poly_cube = 3,
0245                 snyder_poly_octahedron = 4, snyder_poly_dodecahedron = 5,
0246                 snyder_poly_icosahedron = 6
0247             };
0248 
0249             template <typename T>
0250             struct snyder_constants {
0251                 T          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
0252             };
0253 
0254             template <typename T>
0255             inline const snyder_constants<T> * constants()
0256             {
0257                 /* TODO put these in radians to avoid a later conversion */
0258                 static snyder_constants<T> result[] = {
0259                     {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
0260                     {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
0261                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
0262                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
0263                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
0264                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
0265                     {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0}
0266                 };
0267                 return result;
0268             }
0269 
0270             template <typename T>
0271             inline const isea_geo<T> * vertex()
0272             {
0273                 static isea_geo<T> result[] = {
0274                     { 0.0,              deg90_rad<T>()},
0275                     { deg180_rad<T>(),  v_lat},
0276                     {-deg108_rad<T>(),  v_lat},
0277                     {-deg36_rad<T>(),   v_lat},
0278                     { deg36_rad<T>(),   v_lat},
0279                     { deg108_rad<T>(),  v_lat},
0280                     {-deg144_rad<T>(), -v_lat},
0281                     {-deg72_rad<T>(),  -v_lat},
0282                     { 0.0,             -v_lat},
0283                     { deg72_rad<T>(),  -v_lat},
0284                     { deg144_rad<T>(), -v_lat},
0285                     { 0.0,             -deg90_rad<T>()}
0286                 };
0287                 return result;
0288             }
0289 
0290             /* TODO make an isea_pt array of the vertices as well */
0291 
0292             static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
0293 
0294             /* triangle Centers */
0295             template <typename T>
0296             inline const isea_geo<T> * icostriangles()
0297             {
0298                 static isea_geo<T> result[] = {
0299                     { 0.0,              0.0},
0300                     {-deg144_rad<T>(),  e_rad},
0301                     {-deg72_rad<T>(),   e_rad},
0302                     { 0.0,              e_rad},
0303                     { deg72_rad<T>(),   e_rad},
0304                     { deg144_rad<T>(),  e_rad},
0305                     {-deg144_rad<T>(),  f_rad},
0306                     {-deg72_rad<T>(),   f_rad},
0307                     { 0.0,              f_rad},
0308                     { deg72_rad<T>(),   f_rad},
0309                     { deg144_rad<T>(),  f_rad},
0310                     {-deg108_rad<T>(), -f_rad},
0311                     {-deg36_rad<T>(),  -f_rad},
0312                     { deg36_rad<T>(),  -f_rad},
0313                     { deg108_rad<T>(), -f_rad},
0314                     { deg180_rad<T>(), -f_rad},
0315                     {-deg108_rad<T>(), -e_rad},
0316                     {-deg36_rad<T>(),  -e_rad},
0317                     { deg36_rad<T>(),  -e_rad},
0318                     { deg108_rad<T>(), -e_rad},
0319                     { deg180_rad<T>(), -e_rad},
0320                 };
0321                 return result;
0322             }
0323 
0324             template <typename T>
0325             inline T az_adjustment(int triangle)
0326             {
0327                 T          adj;
0328 
0329                 isea_geo<T> v;
0330                 isea_geo<T> c;
0331 
0332                 v = vertex<T>()[tri_v1[triangle]];
0333                 c = icostriangles<T>()[triangle];
0334 
0335                 /* TODO looks like the adjustment is always either 0 or 180 */
0336                 /* at least if you pick your vertex carefully */
0337                 adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
0338                         cos(c.lat) * sin(v.lat)
0339                         - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
0340                 return adj;
0341             }
0342 
0343             template <typename T>
0344             inline isea_pt<T> isea_triangle_xy(int triangle)
0345             {
0346                 isea_pt<T>  c;
0347                 T Rprime = 0.91038328153090290025;
0348 
0349                 triangle = (triangle - 1) % 20;
0350 
0351                 c.x = table_g * ((triangle % 5) - 2) * 2.0;
0352                 if (triangle > 9) {
0353                     c.x += table_g;
0354                 }
0355                 switch (triangle / 5) {
0356                 case 0:
0357                     c.y = 5.0 * table_h;
0358                     break;
0359                 case 1:
0360                     c.y = table_h;
0361                     break;
0362                 case 2:
0363                     c.y = -table_h;
0364                     break;
0365                 case 3:
0366                     c.y = -5.0 * table_h;
0367                     break;
0368                 default:
0369                     /* should be impossible */
0370                     BOOST_THROW_EXCEPTION( projection_exception() );
0371                 };
0372                 c.x *= Rprime;
0373                 c.y *= Rprime;
0374 
0375                 return c;
0376             }
0377 
0378             /* snyder eq 14 */
0379             template <typename T>
0380             inline T sph_azimuth(T const& f_lon, T const& f_lat, T const& t_lon, T const& t_lat)
0381             {
0382                 T          az;
0383 
0384                 az = atan2(cos(t_lat) * sin(t_lon - f_lon),
0385                        cos(f_lat) * sin(t_lat)
0386                        - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
0387                     );
0388                 return az;
0389             }
0390 
0391             /* coord needs to be in radians */
0392             template <typename T>
0393             inline int isea_snyder_forward(isea_geo<T> * ll, isea_pt<T> * out)
0394             {
0395                 static T const two_pi = detail::two_pi<T>();
0396                 static T const d2r = geometry::math::d2r<T>();
0397 
0398                 int             i;
0399 
0400                 /*
0401                  * spherical distance from center of polygon face to any of its
0402                  * vertexes on the globe
0403                  */
0404                 T          g;
0405 
0406                 /*
0407                  * spherical angle between radius vector to center and adjacent edge
0408                  * of spherical polygon on the globe
0409                  */
0410                 T          G;
0411 
0412                 /*
0413                  * plane angle between radius vector to center and adjacent edge of
0414                  * plane polygon
0415                  */
0416                 T          theta;
0417 
0418                 /* additional variables from snyder */
0419                 T          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
0420                                 x, y;
0421 
0422                 /* variables used to store intermediate results */
0423                 T          cot_theta, tan_g, az_offset;
0424 
0425                 /* how many multiples of 60 degrees we adjust the azimuth */
0426                 int             Az_adjust_multiples;
0427 
0428                 snyder_constants<T> c;
0429 
0430                 /*
0431                  * TODO by locality of reference, start by trying the same triangle
0432                  * as last time
0433                  */
0434 
0435                 /* TODO put these constants in as radians to begin with */
0436                 c = constants<T>()[snyder_poly_icosahedron];
0437                 theta = c.theta * d2r;
0438                 g = c.g * d2r;
0439                 G = c.G * d2r;
0440 
0441                 for (i = 1; i <= 20; i++) {
0442                     T          z;
0443                     isea_geo<T> center;
0444 
0445                     center = icostriangles<T>()[i];
0446 
0447                     /* step 1 */
0448                     z = acos(sin(center.lat) * sin(ll->lat)
0449                          + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
0450 
0451                     /* not on this triangle */
0452                     if (z > g + 0.000005) { /* TODO DBL_EPSILON */
0453                         continue;
0454                     }
0455 
0456                     Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);
0457 
0458                     /* step 2 */
0459 
0460                     /* This calculates "some" vertex coordinate */
0461                     az_offset = az_adjustment<T>(i);
0462 
0463                     Az -= az_offset;
0464 
0465                     /* TODO I don't know why we do this.  It's not in snyder */
0466                     /* maybe because we should have picked a better vertex */
0467                     if (Az < 0.0) {
0468                         Az += two_pi;
0469                     }
0470                     /*
0471                      * adjust Az for the point to fall within the range of 0 to
0472                      * 2(90 - theta) or 60 degrees for the hexagon, by
0473                      * and therefore 120 degrees for the triangle
0474                      * of the icosahedron
0475                      * subtracting or adding multiples of 60 degrees to Az and
0476                      * recording the amount of adjustment
0477                      */
0478 
0479                     Az_adjust_multiples = 0;
0480                     while (Az < 0.0) {
0481                         Az += deg120_rad<T>();
0482                         Az_adjust_multiples--;
0483                     }
0484                     while (Az > deg120_rad<T>() + epsilon) {
0485                         Az -= deg120_rad<T>();
0486                         Az_adjust_multiples++;
0487                     }
0488 
0489                     /* step 3 */
0490                     cot_theta = 1.0 / tan(theta);
0491                     tan_g = tan(g);    /* TODO this is a constant */
0492 
0493                     /* Calculate q from eq 9. */
0494                     /* TODO cot_theta is cot(30) */
0495                     q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
0496 
0497                     /* not in this triangle */
0498                     if (z > q + 0.000005) {
0499                         continue;
0500                     }
0501                     /* step 4 */
0502 
0503                     /* Apply equations 5-8 and 10-12 in order */
0504 
0505                     /* eq 5 */
0506                     /* Rprime = 0.9449322893 * R; */
0507                     /* R' in the paper is for the truncated */
0508                     Rprime = 0.91038328153090290025;
0509 
0510                     /* eq 6 */
0511                     H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
0512 
0513                     /* eq 7 */
0514                     /* Ag = (Az + G + H - deg180_rad) * M_PI * R * R / deg180_rad; */
0515                     Ag = Az + G + H - deg180_rad<T>();
0516 
0517                     /* eq 8 */
0518                     Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
0519 
0520                     /* eq 10 */
0521                     /* cot(theta) = 1.73205080756887729355 */
0522                     dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
0523 
0524                     /* eq 11 */
0525                     f = dprime / (2.0 * Rprime * sin(q / 2.0));
0526 
0527                     /* eq 12 */
0528                     rho = 2.0 * Rprime * f * sin(z / 2.0);
0529 
0530                     /*
0531                      * add back the same 60 degree multiple adjustment from step
0532                      * 2 to Azprime
0533                      */
0534 
0535                     Azprime += deg120_rad<T>() * Az_adjust_multiples;
0536 
0537                     /* calculate rectangular coordinates */
0538 
0539                     x = rho * sin(Azprime);
0540                     y = rho * cos(Azprime);
0541 
0542                     /*
0543                      * TODO
0544                      * translate coordinates to the origin for the particular
0545                      * hexagon on the flattened polyhedral map plot
0546                      */
0547 
0548                     out->x = x;
0549                     out->y = y;
0550 
0551                     return i;
0552                 }
0553 
0554                 /*
0555                  * should be impossible, this implies that the coordinate is not on
0556                  * any triangle
0557                  */
0558 
0559                 //fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
0560                 //    ll->lon * geometry::math::r2d<double>(), ll->lat * geometry::math::r2d<double>());
0561                 std::stringstream ss;
0562                 ss << "impossible transform: " << ll->lon * geometry::math::r2d<T>()
0563                    << " " << ll->lat * geometry::math::r2d<T>() << " is not on any triangle.";
0564 
0565                 BOOST_THROW_EXCEPTION( projection_exception(ss.str()) );
0566 
0567                 /* not reached */
0568                 return 0;        /* supresses a warning */
0569             }
0570 
0571             /*
0572              * return the new coordinates of any point in orginal coordinate system.
0573              * Define a point (newNPold) in orginal coordinate system as the North Pole in
0574              * new coordinate system, and the great circle connect the original and new
0575              * North Pole as the lon0 longitude in new coordinate system, given any point
0576              * in orginal coordinate system, this function return the new coordinates.
0577              */
0578 
0579 
0580             /* formula from Snyder, Map Projections: A working manual, p31 */
0581             /*
0582              * old north pole at np in new coordinates
0583              * could be simplified a bit with fewer intermediates
0584              *
0585              * TODO take a result pointer
0586              */
0587             template <typename T>
0588             inline isea_geo<T> snyder_ctran(isea_geo<T> * np, isea_geo<T> * pt)
0589             {
0590                 static T const pi = detail::pi<T>();
0591                 static T const two_pi = detail::two_pi<T>();
0592 
0593                 isea_geo<T> npt;
0594                 T           alpha, phi, lambda, lambda0, beta, lambdap, phip;
0595                 T           sin_phip;
0596                 T           lp_b;    /* lambda prime minus beta */
0597                 T           cos_p, sin_a;
0598 
0599                 phi = pt->lat;
0600                 lambda = pt->lon;
0601                 alpha = np->lat;
0602                 beta = np->lon;
0603                 lambda0 = beta;
0604 
0605                 cos_p = cos(phi);
0606                 sin_a = sin(alpha);
0607 
0608                 /* mpawm 5-7 */
0609                 sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
0610 
0611                 /* mpawm 5-8b */
0612 
0613                 /* use the two argument form so we end up in the right quadrant */
0614                 lp_b = atan2(cos_p * sin(lambda - lambda0),
0615                    (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
0616 
0617                 lambdap = lp_b + beta;
0618 
0619                 /* normalize longitude */
0620                 /* TODO can we just do a modulus ? */
0621                 lambdap = fmod(lambdap, two_pi);
0622                 while (lambdap > pi)
0623                     lambdap -= two_pi;
0624                 while (lambdap < -pi)
0625                     lambdap += two_pi;
0626 
0627                 phip = asin(sin_phip);
0628 
0629                 npt.lat = phip;
0630                 npt.lon = lambdap;
0631 
0632                 return npt;
0633             }
0634 
0635             template <typename T>
0636             inline isea_geo<T> isea_ctran(isea_geo<T> * np, isea_geo<T> * pt, T const& lon0)
0637             {
0638                 static T const pi = detail::pi<T>();
0639                 static T const two_pi = detail::two_pi<T>();
0640 
0641                 isea_geo<T> npt;
0642 
0643                 np->lon += pi;
0644                 npt = snyder_ctran(np, pt);
0645                 np->lon -= pi;
0646 
0647                 npt.lon -= (pi - lon0 + np->lon);
0648 
0649                 /*
0650                  * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
0651                  * vertex 1 these are 180 degrees apart
0652                  */
0653                 npt.lon += pi;
0654                 /* normalize longitude */
0655                 npt.lon = fmod(npt.lon, two_pi);
0656                 while (npt.lon > pi)
0657                     npt.lon -= two_pi;
0658                 while (npt.lon < -pi)
0659                     npt.lon += two_pi;
0660 
0661                 return npt;
0662             }
0663 
0664             /* in radians */
0665 
0666             /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
0667 
0668             template <typename T>
0669             inline int isea_grid_init(isea_dgg<T> * g)
0670             {
0671                 if (!g)
0672                     return 0;
0673 
0674                 //g->polyhedron = isea_icosahedron;
0675                 g->o_lat = isea_std_lat;
0676                 g->o_lon = isea_std_lon;
0677                 g->o_az = 0.0;
0678                 g->aperture = 4;
0679                 g->resolution = 6;
0680                 g->radius = 1.0;
0681                 //g->topology = isea_hexagon;
0682 
0683                 return 1;
0684             }
0685 
0686             template <typename T>
0687             inline int isea_orient_isea(isea_dgg<T> * g)
0688             {
0689                 if (!g)
0690                     return 0;
0691                 g->o_lat = isea_std_lat;
0692                 g->o_lon = isea_std_lon;
0693                 g->o_az = 0.0;
0694                 return 1;
0695             }
0696 
0697             template <typename T>
0698             inline int isea_orient_pole(isea_dgg<T> * g)
0699             {
0700                 static T const half_pi = detail::half_pi<T>();
0701 
0702                 if (!g)
0703                     return 0;
0704                 g->o_lat = half_pi;
0705                 g->o_lon = 0.0;
0706                 g->o_az = 0;
0707                 return 1;
0708             }
0709 
0710             template <typename T>
0711             inline int isea_transform(isea_dgg<T> * g, isea_geo<T> * in,
0712                                       isea_pt<T> * out)
0713             {
0714                 isea_geo<T> i, pole;
0715                 int         tri;
0716 
0717                 pole.lat = g->o_lat;
0718                 pole.lon = g->o_lon;
0719 
0720                 i = isea_ctran(&pole, in, g->o_az);
0721 
0722                 tri = isea_snyder_forward(&i, out);
0723                 out->x *= g->radius;
0724                 out->y *= g->radius;
0725                 g->triangle = tri;
0726 
0727                 return tri;
0728             }
0729 
0730 
0731             template <typename T>
0732             inline void isea_rotate(isea_pt<T> * pt, T const& degrees)
0733             {
0734                 static T const d2r = geometry::math::d2r<T>();
0735                 static T const two_pi = detail::two_pi<T>();
0736 
0737                 T          rad;
0738 
0739                 T          x, y;
0740 
0741                 rad = -degrees * d2r;
0742                 while (rad >= two_pi) rad -= two_pi;
0743                 while (rad <= -two_pi) rad += two_pi;
0744 
0745                 x = pt->x * cos(rad) + pt->y * sin(rad);
0746                 y = -pt->x * sin(rad) + pt->y * cos(rad);
0747 
0748                 pt->x = x;
0749                 pt->y = y;
0750             }
0751 
0752             template <typename T>
0753             inline int isea_tri_plane(int tri, isea_pt<T> *pt, T const& radius)
0754             {
0755                 isea_pt<T> tc; /* center of triangle */
0756 
0757                 if (downtri(tri)) {
0758                     isea_rotate(pt, 180.0);
0759                 }
0760                 tc = isea_triangle_xy<T>(tri);
0761                 tc.x *= radius;
0762                 tc.y *= radius;
0763                 pt->x += tc.x;
0764                 pt->y += tc.y;
0765 
0766                 return tri;
0767             }
0768 
0769             /* convert projected triangle coords to quad xy coords, return quad number */
0770             template <typename T>
0771             inline int isea_ptdd(int tri, isea_pt<T> *pt)
0772             {
0773                 int             downtri, quad;
0774 
0775                 downtri = (((tri - 1) / 5) % 2 == 1);
0776                 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
0777 
0778                 isea_rotate(pt, downtri ? 240.0 : 60.0);
0779                 if (downtri) {
0780                     pt->x += 0.5;
0781                     /* pt->y += cos(30.0 * M_PI / 180.0); */
0782                     pt->y += .86602540378443864672;
0783                 }
0784                 return quad;
0785             }
0786 
0787             template <typename T>
0788             inline int isea_dddi_ap3odd(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di)
0789             {
0790                 static T const pi = detail::pi<T>();
0791 
0792                 isea_pt<T> v;
0793                 T          hexwidth;
0794                 T          sidelength;    /* in hexes */
0795                 int        d, i;
0796                 int        maxcoord;
0797                 hex        h;
0798 
0799                 /* This is the number of hexes from apex to base of a triangle */
0800                 sidelength = (math::pow(T(2), g->resolution) + T(1)) / T(2);
0801 
0802                 /* apex to base is cos(30deg) */
0803                 hexwidth = cos(pi / 6.0) / sidelength;
0804 
0805                 /* TODO I think sidelength is always x.5, so
0806                  * (int)sidelength * 2 + 1 might be just as good
0807                  */
0808                 maxcoord = (int) (sidelength * 2.0 + 0.5);
0809 
0810                 v = *pt;
0811                 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
0812                 h.iso = 0;
0813                 hex_iso(&h);
0814 
0815                 d = h.x - h.z;
0816                 i = h.x + h.y + h.y;
0817 
0818                 /*
0819                  * you want to test for max coords for the next quad in the same
0820                  * "row" first to get the case where both are max
0821                  */
0822                 if (quad <= 5) {
0823                     if (d == 0 && i == maxcoord) {
0824                         /* north pole */
0825                         quad = 0;
0826                         d = 0;
0827                         i = 0;
0828                     } else if (i == maxcoord) {
0829                         /* upper right in next quad */
0830                         quad += 1;
0831                         if (quad == 6)
0832                             quad = 1;
0833                         i = maxcoord - d;
0834                         d = 0;
0835                     } else if (d == maxcoord) {
0836                         /* lower right in quad to lower right */
0837                         quad += 5;
0838                         d = 0;
0839                     }
0840                 } else if (quad >= 6) {
0841                     if (i == 0 && d == maxcoord) {
0842                         /* south pole */
0843                         quad = 11;
0844                         d = 0;
0845                         i = 0;
0846                     } else if (d == maxcoord) {
0847                         /* lower right in next quad */
0848                         quad += 1;
0849                         if (quad == 11)
0850                             quad = 6;
0851                         d = maxcoord - i;
0852                         i = 0;
0853                     } else if (i == maxcoord) {
0854                         /* upper right in quad to upper right */
0855                         quad = (quad - 4) % 5;
0856                         i = 0;
0857                     }
0858                 }
0859 
0860                 di->x = d;
0861                 di->y = i;
0862 
0863                 g->quad = quad;
0864                 return quad;
0865             }
0866 
0867             template <typename T>
0868             inline int isea_dddi(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di)
0869             {
0870                 isea_pt<T> v;
0871                 T          hexwidth;
0872                 int        sidelength;    /* in hexes */
0873                 hex        h;
0874 
0875                 if (g->aperture == 3 && g->resolution % 2 != 0) {
0876                     return isea_dddi_ap3odd(g, quad, pt, di);
0877                 }
0878                 /* todo might want to do this as an iterated loop */
0879                 if (g->aperture >0) {
0880                     sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5));
0881                 } else {
0882                     sidelength = g->resolution;
0883                 }
0884 
0885                 hexwidth = 1.0 / sidelength;
0886 
0887                 v = *pt;
0888                 isea_rotate(&v, -30.0);
0889                 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
0890                 h.iso = 0;
0891                 hex_iso(&h);
0892 
0893                 /* we may actually be on another quad */
0894                 if (quad <= 5) {
0895                     if (h.x == 0 && h.z == -sidelength) {
0896                         /* north pole */
0897                         quad = 0;
0898                         h.z = 0;
0899                         h.y = 0;
0900                         h.x = 0;
0901                     } else if (h.z == -sidelength) {
0902                         quad = quad + 1;
0903                         if (quad == 6)
0904                             quad = 1;
0905                         h.y = sidelength - h.x;
0906                         h.z = h.x - sidelength;
0907                         h.x = 0;
0908                     } else if (h.x == sidelength) {
0909                         quad += 5;
0910                         h.y = -h.z;
0911                         h.x = 0;
0912                     }
0913                 } else if (quad >= 6) {
0914                     if (h.z == 0 && h.x == sidelength) {
0915                         /* south pole */
0916                         quad = 11;
0917                         h.x = 0;
0918                         h.y = 0;
0919                         h.z = 0;
0920                     } else if (h.x == sidelength) {
0921                         quad = quad + 1;
0922                         if (quad == 11)
0923                             quad = 6;
0924                         h.x = h.y + sidelength;
0925                         h.y = 0;
0926                         h.z = -h.x;
0927                     } else if (h.y == -sidelength) {
0928                         quad -= 4;
0929                         h.y = 0;
0930                         h.z = -h.x;
0931                     }
0932                 }
0933                 di->x = h.x;
0934                 di->y = -h.z;
0935 
0936                 g->quad = quad;
0937                 return quad;
0938             }
0939 
0940             template <typename T>
0941             inline int isea_ptdi(isea_dgg<T> *g, int tri, isea_pt<T> *pt,
0942                                  isea_pt<T> *di)
0943             {
0944                 isea_pt<T> v;
0945                 int        quad;
0946 
0947                 v = *pt;
0948                 quad = isea_ptdd(tri, &v);
0949                 quad = isea_dddi(g, quad, &v, di);
0950                 return quad;
0951             }
0952 
0953             /* q2di to seqnum */
0954             template <typename T>
0955             inline int isea_disn(isea_dgg<T> *g, int quad, isea_pt<T> *di)
0956             {
0957                 int             sidelength;
0958                 int             sn, height;
0959                 int             hexes;
0960 
0961                 if (quad == 0) {
0962                     g->serial = 1;
0963                     return g->serial;
0964                 }
0965                 /* hexes in a quad */
0966                 hexes = (int) (math::pow(T(g->aperture), T(g->resolution)) + T(0.5));
0967                 if (quad == 11) {
0968                     g->serial = 1 + 10 * hexes + 1;
0969                     return g->serial;
0970                 }
0971                 if (g->aperture == 3 && g->resolution % 2 == 1) {
0972                     height = (int) (math::pow(T(g->aperture), T((g->resolution - 1) / T(2))));
0973                     sn = ((int) di->x) * height;
0974                     sn += ((int) di->y) / height;
0975                     sn += (quad - 1) * hexes;
0976                     sn += 2;
0977                 } else {
0978                     sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5));
0979                     sn = (int) ((quad - 1) * hexes + sidelength * di->x + di->y + 2);
0980                 }
0981 
0982                 g->serial = sn;
0983                 return sn;
0984             }
0985 
0986             /* TODO just encode the quad in the d or i coordinate
0987              * quad is 0-11, which can be four bits.
0988              * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
0989              */
0990             /* convert a q2di to global hex coord */
0991             template <typename T>
0992             inline int isea_hex(isea_dgg<T> *g, int tri, isea_pt<T> *pt,
0993                                 isea_pt<T> *hex)
0994             {
0995                 isea_pt<T> v;
0996 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME
0997                 int sidelength;
0998                 int d, i, x, y;
0999 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME
1000                 int quad;
1001 
1002                 quad = isea_ptdi(g, tri, pt, &v);
1003 
1004                 hex->x = ((int)v.x << 4) + quad;
1005                 hex->y = v.y;
1006 
1007                 return 1;
1008 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME
1009                 d = (int)v.x;
1010                 i = (int)v.y;
1011 
1012                 /* Aperture 3 odd resolutions */
1013                 if (g->aperture == 3 && g->resolution % 2 != 0) {
1014                     int offset = (int)(pow(T(3.0), T(g->resolution - 1)) + 0.5);
1015 
1016                     d += offset * ((g->quad-1) % 5);
1017                     i += offset * ((g->quad-1) % 5);
1018 
1019                     if (quad == 0) {
1020                         d = 0;
1021                         i = offset;
1022                     } else if (quad == 11) {
1023                         d = 2 * offset;
1024                         i = 0;
1025                     } else if (quad > 5) {
1026                         d += offset;
1027                     }
1028 
1029                     x = (2*d - i) /3;
1030                     y = (2*i - d) /3;
1031 
1032                     hex->x = x + offset / 3;
1033                     hex->y = y + 2 * offset / 3;
1034                     return 1;
1035                 }
1036 
1037                 /* aperture 3 even resolutions and aperture 4 */
1038                 sidelength = (int) (pow(T(g->aperture), T(g->resolution / 2.0)) + 0.5);
1039                 if (g->quad == 0) {
1040                     hex->x = 0;
1041                     hex->y = sidelength;
1042                 } else if (g->quad == 11) {
1043                     hex->x = sidelength * 2;
1044                     hex->y = 0;
1045                 } else {
1046                     hex->x = d + sidelength * ((g->quad-1) % 5);
1047                     if (g->quad > 5) hex->x += sidelength;
1048                     hex->y = i + sidelength * ((g->quad-1) % 5);
1049                 }
1050 
1051                 return 1;
1052 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME
1053             }
1054 
1055             template <typename T>
1056             inline isea_pt<T> isea_forward(isea_dgg<T> *g, isea_geo<T> *in)
1057             {
1058                 int        tri;
1059                 isea_pt<T> out, coord;
1060 
1061                 tri = isea_transform(g, in, &out);
1062 
1063                 if (g->output == isea_addr_plane) {
1064                     isea_tri_plane(tri, &out, g->radius);
1065                     return out;
1066                 }
1067 
1068                 /* convert to isea standard triangle size */
1069                 out.x = out.x / g->radius * isea_scale;
1070                 out.y = out.y / g->radius * isea_scale;
1071                 out.x += 0.5;
1072                 out.y += 2.0 * .14433756729740644112;
1073 
1074                 switch (g->output) {
1075                     case isea_addr_projtri:
1076                         /* nothing to do, already in projected triangle */
1077                         break;
1078                     case isea_addr_vertex2dd:
1079                         g->quad = isea_ptdd(tri, &out);
1080                         break;
1081                     case isea_addr_q2dd:
1082                         /* Same as above, we just don't print as much */
1083                         g->quad = isea_ptdd(tri, &out);
1084                         break;
1085                     case isea_addr_q2di:
1086                         g->quad = isea_ptdi(g, tri, &out, &coord);
1087                         return coord;
1088                         break;
1089                     case isea_addr_seqnum:
1090                         isea_ptdi(g, tri, &out, &coord);
1091                         /* disn will set g->serial */
1092                         isea_disn(g, g->quad, &coord);
1093                         return coord;
1094                         break;
1095                     case isea_addr_hex:
1096                         isea_hex(g, tri, &out, &coord);
1097                         return coord;
1098                         break;
1099                     default:
1100                         // isea_addr_plane handled above
1101                         BOOST_GEOMETRY_ASSERT(false);
1102                         break;
1103                 }
1104 
1105                 return out;
1106             }
1107             /*
1108              * Proj 4 integration code follows
1109              */
1110 
1111             template <typename T>
1112             struct par_isea
1113             {
1114                 isea_dgg<T> dgg;
1115             };
1116 
1117             template <typename T, typename Parameters>
1118             struct base_isea_spheroid
1119             {
1120                 par_isea<T> m_proj_parm;
1121 
1122                 // FORWARD(s_forward)
1123                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
1124                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
1125                 {
1126                     isea_pt<T> out;
1127                     isea_geo<T> in;
1128 
1129                     in.lon = lp_lon;
1130                     in.lat = lp_lat;
1131 
1132                     isea_dgg<T> copy = this->m_proj_parm.dgg;
1133                     out = isea_forward(&copy, &in);
1134 
1135                     xy_x = out.x;
1136                     xy_y = out.y;
1137                 }
1138 
1139                 static inline std::string get_name()
1140                 {
1141                     return "isea_spheroid";
1142                 }
1143 
1144             };
1145 
1146             template <typename T>
1147             inline void isea_orient_init(srs::detail::proj4_parameters const& params,
1148                                          par_isea<T>& proj_parm)
1149             {
1150                 std::string opt = pj_get_param_s(params, "orient");
1151                 if (! opt.empty()) {
1152                     if (opt == std::string("isea")) {
1153                         isea_orient_isea(&proj_parm.dgg);
1154                     } else if (opt == std::string("pole")) {
1155                         isea_orient_pole(&proj_parm.dgg);
1156                     } else {
1157                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1158                     }
1159                 }
1160             }
1161 
1162             template <typename T>
1163             inline void isea_orient_init(srs::dpar::parameters<T> const& params,
1164                                          par_isea<T>& proj_parm)
1165             {
1166                 typename srs::dpar::parameters<T>::const_iterator
1167                     it = pj_param_find(params, srs::dpar::orient);
1168                 if (it != params.end()) {
1169                     srs::dpar::value_orient o = static_cast<srs::dpar::value_orient>(it->template get_value<int>());
1170                     if (o == srs::dpar::orient_isea) {
1171                         isea_orient_isea(&proj_parm.dgg);
1172                     } else if (o == srs::dpar::orient_pole) {
1173                         isea_orient_pole(&proj_parm.dgg);
1174                     } else {
1175                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1176                     }
1177                 }
1178             }
1179 
1180             template <typename T>
1181             inline void isea_mode_init(srs::detail::proj4_parameters const& params,
1182                                        par_isea<T>& proj_parm)
1183             {
1184                 std::string opt = pj_get_param_s(params, "mode");
1185                 if (! opt.empty()) {
1186                     if (opt == std::string("plane")) {
1187                         proj_parm.dgg.output = isea_addr_plane;
1188                     } else if (opt == std::string("di")) {
1189                         proj_parm.dgg.output = isea_addr_q2di;
1190                     } else if (opt == std::string("dd")) {
1191                         proj_parm.dgg.output = isea_addr_q2dd;
1192                     } else if (opt == std::string("hex")) {
1193                         proj_parm.dgg.output = isea_addr_hex;
1194                     } else {
1195                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1196                     }
1197                 }
1198             }
1199 
1200             template <typename T>
1201             inline void isea_mode_init(srs::dpar::parameters<T> const& params,
1202                                        par_isea<T>& proj_parm)
1203             {
1204                 typename srs::dpar::parameters<T>::const_iterator
1205                     it = pj_param_find(params, srs::dpar::mode);
1206                 if (it != params.end()) {
1207                     srs::dpar::value_mode m = static_cast<srs::dpar::value_mode>(it->template get_value<int>());
1208                     if (m == srs::dpar::mode_plane) {
1209                         proj_parm.dgg.output = isea_addr_plane;
1210                     } else if (m == srs::dpar::mode_di) {
1211                         proj_parm.dgg.output = isea_addr_q2di;
1212                     } else if (m == srs::dpar::mode_dd) {
1213                         proj_parm.dgg.output = isea_addr_q2dd;
1214                     } else if (m == srs::dpar::mode_hex) {
1215                         proj_parm.dgg.output = isea_addr_hex;
1216                     } else {
1217                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1218                     }
1219                 }
1220             }
1221 
1222             // Icosahedral Snyder Equal Area
1223             template <typename Params, typename T>
1224             inline void setup_isea(Params const& params, par_isea<T>& proj_parm)
1225             {
1226                 std::string opt;
1227 
1228                 isea_grid_init(&proj_parm.dgg);
1229 
1230                 proj_parm.dgg.output = isea_addr_plane;
1231             /*        proj_parm.dgg.radius = par.a; / * otherwise defaults to 1 */
1232                 /* calling library will scale, I think */
1233 
1234                 isea_orient_init(params, proj_parm);
1235 
1236                 pj_param_r<srs::spar::azi>(params, "azi", srs::dpar::azi, proj_parm.dgg.o_az);
1237                 pj_param_r<srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0, proj_parm.dgg.o_lon);
1238                 pj_param_r<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0, proj_parm.dgg.o_lat);
1239                 // TODO: this parameter is set below second time
1240                 pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture);
1241                 // TODO: this parameter is set below second time
1242                 pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution);
1243 
1244                 isea_mode_init(params, proj_parm);
1245 
1246                 // TODO: pj_param_exists -> pj_get_param_b ?
1247                 if (pj_param_exists<srs::spar::rescale>(params, "rescale", srs::dpar::rescale)) {
1248                     proj_parm.dgg.radius = isea_scale;
1249                 }
1250 
1251                 if (pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution)) {
1252                     /* empty */
1253                 } else {
1254                     proj_parm.dgg.resolution = 4;
1255                 }
1256 
1257                 if (pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture)) {
1258                     /* empty */
1259                 } else {
1260                     proj_parm.dgg.aperture = 3;
1261                 }
1262             }
1263 
1264     }} // namespace detail::isea
1265     #endif // doxygen
1266 
1267     /*!
1268         \brief Icosahedral Snyder Equal Area projection
1269         \ingroup projections
1270         \tparam Geographic latlong point type
1271         \tparam Cartesian xy point type
1272         \tparam Parameters parameter type
1273         \par Projection characteristics
1274          - Spheroid
1275         \par Projection parameters
1276          - orient (string)
1277          - azi: Azimuth (or Gamma) (degrees)
1278          - lon_0: Central meridian (degrees)
1279          - lat_0: Latitude of origin (degrees)
1280          - aperture (integer)
1281          - resolution (integer)
1282          - mode (string)
1283          - rescale
1284         \par Example
1285         \image html ex_isea.gif
1286     */
1287     template <typename T, typename Parameters>
1288     struct isea_spheroid : public detail::isea::base_isea_spheroid<T, Parameters>
1289     {
1290         template <typename Params>
1291         inline isea_spheroid(Params const& params, Parameters const& )
1292         {
1293             detail::isea::setup_isea(params, this->m_proj_parm);
1294         }
1295     };
1296 
1297     #ifndef DOXYGEN_NO_DETAIL
1298     namespace detail
1299     {
1300 
1301         // Static projection
1302         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_F(srs::spar::proj_isea, isea_spheroid)
1303 
1304         // Factory entry(s)
1305         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_F(isea_entry, isea_spheroid)
1306 
1307         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(isea_init)
1308         {
1309             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(isea, isea_entry)
1310         }
1311 
1312     } // namespace detail
1313     #endif // doxygen
1314 
1315 } // namespace projections
1316 
1317 }} // namespace boost::geometry
1318 
1319 #endif // BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
1320
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!