Back to home page

EIC code displayed by LXR
 
 

    

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

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 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028 
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031 
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039 
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
0042 
0043 #include <boost/config.hpp>
0044 #include <boost/geometry/util/math.hpp>
0045 #include <boost/math/special_functions/hypot.hpp>
0046 
0047 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0048 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0049 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0050 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0051 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
0052 #include <boost/geometry/srs/projections/impl/projects.hpp>
0053 
0054 namespace boost { namespace geometry
0055 {
0056 
0057 namespace projections
0058 {
0059     #ifndef DOXYGEN_NO_DETAIL
0060     namespace detail { namespace stere
0061     {
0062             static const double epsilon10 = 1.e-10;
0063             static const double tolerance = 1.e-8;
0064             static const int n_iter = 8;
0065             static const double conv_tolerance = 1.e-10;
0066 
0067             enum mode_type {
0068                 s_pole = 0,
0069                 n_pole = 1,
0070                 obliq  = 2,
0071                 equit  = 3
0072             };
0073 
0074             template <typename T>
0075             struct par_stere
0076             {
0077                 T   phits;
0078                 T   sinX1;
0079                 T   cosX1;
0080                 T   akm1;
0081                 mode_type mode;
0082                 bool variant_c;
0083             };
0084 
0085             template <typename T>
0086             inline T ssfn_(T const& phit, T sinphi, T const& eccen)
0087             {
0088                 static const T half_pi = detail::half_pi<T>();
0089 
0090                 sinphi *= eccen;
0091                 return (tan (.5 * (half_pi + phit)) *
0092                    math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * eccen));
0093             }
0094 
0095             template <typename T, typename Parameters>
0096             struct base_stere_ellipsoid
0097             {
0098                 par_stere<T> m_proj_parm;
0099 
0100                 // FORWARD(e_forward)  ellipsoid
0101                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0102                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0103                 {
0104                     static const T half_pi = detail::half_pi<T>();
0105 
0106                     T coslam, sinlam, sinX=0.0, cosX=0.0, X, A = 0.0, sinphi;
0107 
0108                     coslam = cos(lp_lon);
0109                     sinlam = sin(lp_lon);
0110                     sinphi = sin(lp_lat);
0111                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
0112                         sinX = sin(X = 2. * atan(ssfn_(lp_lat, sinphi, par.e)) - half_pi);
0113                         cosX = cos(X);
0114                     }
0115                     switch (this->m_proj_parm.mode) {
0116                     case obliq:
0117                         A = this->m_proj_parm.akm1 / (this->m_proj_parm.cosX1 * (1. + this->m_proj_parm.sinX1 * sinX +
0118                            this->m_proj_parm.cosX1 * cosX * coslam));
0119                         xy_y = A * (this->m_proj_parm.cosX1 * sinX - this->m_proj_parm.sinX1 * cosX * coslam);
0120                         goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */
0121 
0122                     case equit:
0123                         // TODO: calculate denominator once
0124                         /* avoid zero division */
0125                         if (1. + cosX * coslam == 0.0) {
0126                             xy_y = HUGE_VAL;
0127                         } else {
0128                             A = this->m_proj_parm.akm1 / (1. + cosX * coslam);
0129                             xy_y = A * sinX;
0130                         }
0131                 xmul:
0132                         xy_x = A * cosX;
0133                         break;
0134 
0135                     case s_pole:
0136                         lp_lat = -lp_lat;
0137                         coslam = - coslam;
0138                         sinphi = -sinphi;
0139                         BOOST_FALLTHROUGH;
0140                     case n_pole:
0141                         // see IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2
0142                         // December 2021 pg. 82
0143                         if( m_proj_parm.variant_c )
0144                         {
0145                             auto t = pj_tsfn(lp_lat, sinphi, par.e);
0146                             auto tf = pj_tsfn(this->m_proj_parm.phits,
0147                                                 sin(this->m_proj_parm.phits),
0148                                                 par.e);
0149                             xy_x = this->m_proj_parm.akm1 * t;
0150                             auto mf = this->m_proj_parm.akm1 * tf;
0151                             xy_y = - xy_x * coslam - mf;
0152                         } else {
0153                             xy_x = this->m_proj_parm.akm1 * pj_tsfn(lp_lat, sinphi, par.e);
0154                             xy_y = - xy_x * coslam;
0155                         }
0156                         break;
0157                     }
0158 
0159                     xy_x = xy_x * sinlam;
0160                 }
0161 
0162                 // INVERSE(e_inverse)  ellipsoid
0163                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0164                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0165                 {
0166                     static const T half_pi = detail::half_pi<T>();
0167 
0168                     T cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
0169                     T mf = 0;
0170                     int i;
0171 
0172                     rho = boost::math::hypot(xy_x, xy_y);
0173                     switch (this->m_proj_parm.mode) {
0174                     case obliq:
0175                     case equit:
0176                         cosphi = cos( tp = 2. * atan2(rho * this->m_proj_parm.cosX1 , this->m_proj_parm.akm1) );
0177                         sinphi = sin(tp);
0178                         if( rho == 0.0 )
0179                             phi_l = asin(cosphi * this->m_proj_parm.sinX1);
0180                         else
0181                             phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));
0182 
0183                         tp = tan(.5 * (half_pi + phi_l));
0184                         xy_x *= sinphi;
0185                         xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
0186                         halfpi = half_pi;
0187                         halfe = .5 * par.e;
0188                         break;
0189                     case n_pole:
0190                         xy_y = -xy_y;
0191                         BOOST_FALLTHROUGH;
0192                     case s_pole:
0193                         // see IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2
0194                         // December 2021 pg. 82
0195                         if( m_proj_parm.variant_c )
0196                         {
0197                             auto tf = pj_tsfn(this->m_proj_parm.phits,
0198                                               sin(this->m_proj_parm.phits),
0199                                               par.e);
0200                             mf = this->m_proj_parm.akm1 * tf;
0201                             rho = boost::math::hypot(xy_x, xy_y + mf);
0202                         }
0203                         phi_l = half_pi - 2. * atan(tp = - rho / this->m_proj_parm.akm1);
0204                         halfpi = -half_pi;
0205                         halfe = -.5 * par.e;
0206                         break;
0207                     }
0208                     for (i = n_iter; i--; phi_l = lp_lat) {
0209                         sinphi = par.e * sin(phi_l);
0210                         lp_lat = T(2) * atan(tp * math::pow((T(1)+sinphi)/(T(1)-sinphi), halfe)) - halfpi;
0211                         if (fabs(phi_l - lp_lat) < conv_tolerance) {
0212                             if (this->m_proj_parm.mode == s_pole)
0213                                 lp_lat = -lp_lat;
0214                             lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y + mf);
0215                             return;
0216                         }
0217                     }
0218                     BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0219                 }
0220 
0221                 static inline std::string get_name()
0222                 {
0223                     return "stere_ellipsoid";
0224                 }
0225 
0226             };
0227 
0228             template <typename T, typename Parameters>
0229             struct base_stere_spheroid
0230             {
0231                 par_stere<T> m_proj_parm;
0232 
0233                 // FORWARD(s_forward)  spheroid
0234                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0235                 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0236                 {
0237                     static const T fourth_pi = detail::fourth_pi<T>();
0238                     static const T half_pi = detail::half_pi<T>();
0239 
0240                     T  sinphi, cosphi, coslam, sinlam;
0241 
0242                     sinphi = sin(lp_lat);
0243                     cosphi = cos(lp_lat);
0244                     coslam = cos(lp_lon);
0245                     sinlam = sin(lp_lon);
0246                     switch (this->m_proj_parm.mode) {
0247                     case equit:
0248                         xy_y = 1. + cosphi * coslam;
0249                         goto oblcon;
0250                     case obliq:
0251                         xy_y = 1. + this->m_proj_parm.sinX1 * sinphi + this->m_proj_parm.cosX1 * cosphi * coslam;
0252                 oblcon:
0253                         if (xy_y <= epsilon10) {
0254                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0255                         }
0256                         xy_x = (xy_y = this->m_proj_parm.akm1 / xy_y) * cosphi * sinlam;
0257                         xy_y *= (this->m_proj_parm.mode == equit) ? sinphi :
0258                            this->m_proj_parm.cosX1 * sinphi - this->m_proj_parm.sinX1 * cosphi * coslam;
0259                         break;
0260                     case n_pole:
0261                         coslam = - coslam;
0262                         lp_lat = - lp_lat;
0263                         BOOST_FALLTHROUGH;
0264                     case s_pole:
0265                         if (fabs(lp_lat - half_pi) < tolerance) {
0266                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0267                         }
0268                         xy_x = sinlam * ( xy_y = this->m_proj_parm.akm1 * tan(fourth_pi + .5 * lp_lat) );
0269                         xy_y *= coslam;
0270                         break;
0271                     }
0272                 }
0273 
0274                 // INVERSE(s_inverse)  spheroid
0275                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0276                 inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0277                 {
0278                     T  c, rh, sinc, cosc;
0279 
0280                     sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
0281                     cosc = cos(c);
0282                     lp_lon = 0.;
0283 
0284                     switch (this->m_proj_parm.mode) {
0285                     case equit:
0286                         if (fabs(rh) <= epsilon10)
0287                             lp_lat = 0.;
0288                         else
0289                             lp_lat = asin(xy_y * sinc / rh);
0290                         if (cosc != 0. || xy_x != 0.)
0291                             lp_lon = atan2(xy_x * sinc, cosc * rh);
0292                         break;
0293                     case obliq:
0294                         if (fabs(rh) <= epsilon10)
0295                             lp_lat = par.phi0;
0296                         else
0297                             lp_lat = asin(cosc * this->m_proj_parm.sinX1 + xy_y * sinc * this->m_proj_parm.cosX1 / rh);
0298                         if ((c = cosc - this->m_proj_parm.sinX1 * sin(lp_lat)) != 0. || xy_x != 0.)
0299                             lp_lon = atan2(xy_x * sinc * this->m_proj_parm.cosX1, c * rh);
0300                         break;
0301                     case n_pole:
0302                         xy_y = -xy_y;
0303                         BOOST_FALLTHROUGH;
0304                     case s_pole:
0305                         if (fabs(rh) <= epsilon10)
0306                             lp_lat = par.phi0;
0307                         else
0308                             lp_lat = asin(this->m_proj_parm.mode == s_pole ? - cosc : cosc);
0309                         lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
0310                         break;
0311                     }
0312                 }
0313 
0314                 static inline std::string get_name()
0315                 {
0316                     return "stere_spheroid";
0317                 }
0318 
0319             };
0320 
0321             template <typename Parameters, typename T>
0322             inline void setup(Parameters const& par, par_stere<T>& proj_parm)  /* general initialization */
0323             {
0324                 static const T fourth_pi = detail::fourth_pi<T>();
0325                 static const T half_pi = detail::half_pi<T>();
0326 
0327                 T t;
0328 
0329                 if (fabs((t = fabs(par.phi0)) - half_pi) < epsilon10)
0330                     proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
0331                 else
0332                     proj_parm.mode = t > epsilon10 ? obliq : equit;
0333                 proj_parm.phits = fabs(proj_parm.phits);
0334 
0335                 if (par.es != 0.0) {
0336                     T X;
0337 
0338                     switch (proj_parm.mode) {
0339                     case n_pole:
0340                     case s_pole:
0341                         if (fabs(proj_parm.phits - half_pi) < epsilon10)
0342                             proj_parm.akm1 = 2. * par.k0 /
0343                                sqrt(math::pow(T(1)+par.e,T(1)+par.e)*math::pow(T(1)-par.e,T(1)-par.e));
0344                         else {
0345                             proj_parm.akm1 = cos(proj_parm.phits) /
0346                                pj_tsfn(proj_parm.phits, t = sin(proj_parm.phits), par.e);
0347                             t *= par.e;
0348                             proj_parm.akm1 /= sqrt(1. - t * t);
0349                         }
0350                         break;
0351                     case equit:
0352                     case obliq:
0353                         t = sin(par.phi0);
0354                         X = 2. * atan(ssfn_(par.phi0, t, par.e)) - half_pi;
0355                         t *= par.e;
0356                         proj_parm.akm1 = 2. * par.k0 * cos(par.phi0) / sqrt(1. - t * t);
0357                         proj_parm.sinX1 = sin(X);
0358                         proj_parm.cosX1 = cos(X);
0359                         break;
0360                     }
0361                 } else {
0362                     switch (proj_parm.mode) {
0363                     case obliq:
0364                         proj_parm.sinX1 = sin(par.phi0);
0365                         proj_parm.cosX1 = cos(par.phi0);
0366                         BOOST_FALLTHROUGH;
0367                     case equit:
0368                         proj_parm.akm1 = 2. * par.k0;
0369                         break;
0370                     case s_pole:
0371                     case n_pole:
0372                         proj_parm.akm1 = fabs(proj_parm.phits - half_pi) >= epsilon10 ?
0373                            cos(proj_parm.phits) / tan(fourth_pi - .5 * proj_parm.phits) :
0374                            2. * par.k0 ;
0375                         break;
0376                     }
0377                 }
0378             }
0379 
0380 
0381             // Stereographic
0382             template <typename Params, typename Parameters, typename T>
0383             inline void setup_stere(Params const& params, Parameters const& par, par_stere<T>& proj_parm)
0384             {
0385                 static const T half_pi = detail::half_pi<T>();
0386 
0387                 if (! pj_param_r<srs::spar::lat_ts>(params, "lat_ts", srs::dpar::lat_ts, proj_parm.phits))
0388                     proj_parm.phits = half_pi;
0389 
0390                 proj_parm.variant_c = false;
0391                 if (pj_param_exists<srs::spar::variant_c>(params, "variant_c", srs::dpar::variant_c))
0392                     proj_parm.variant_c = true;
0393 
0394                 setup(par, proj_parm);
0395             }
0396 
0397             // Universal Polar Stereographic
0398             template <typename Params, typename Parameters, typename T>
0399             inline void setup_ups(Params const& params, Parameters& par, par_stere<T>& proj_parm)
0400             {
0401                 static const T half_pi = detail::half_pi<T>();
0402 
0403                 /* International Ellipsoid */
0404                 par.phi0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi: half_pi;
0405                 if (par.es == 0.0) {
0406                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
0407                 }
0408                 par.k0 = .994;
0409                 par.x0 = 2000000.;
0410                 par.y0 = 2000000.;
0411                 proj_parm.phits = half_pi;
0412                 par.lam0 = 0.;
0413 
0414                 setup(par, proj_parm);
0415             }
0416 
0417     }} // namespace detail::stere
0418     #endif // doxygen
0419 
0420     /*!
0421         \brief Stereographic projection
0422         \ingroup projections
0423         \tparam Geographic latlong point type
0424         \tparam Cartesian xy point type
0425         \tparam Parameters parameter type
0426         \par Projection characteristics
0427          - Azimuthal
0428          - Spheroid
0429          - Ellipsoid
0430         \par Projection parameters
0431          - lat_ts: Latitude of true scale (degrees)
0432         \par Example
0433         \image html ex_stere.gif
0434     */
0435     template <typename T, typename Parameters>
0436     struct stere_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
0437     {
0438         template <typename Params>
0439         inline stere_ellipsoid(Params const& params, Parameters const& par)
0440         {
0441             detail::stere::setup_stere(params, par, this->m_proj_parm);
0442         }
0443     };
0444 
0445     /*!
0446         \brief Stereographic projection
0447         \ingroup projections
0448         \tparam Geographic latlong point type
0449         \tparam Cartesian xy point type
0450         \tparam Parameters parameter type
0451         \par Projection characteristics
0452          - Azimuthal
0453          - Spheroid
0454          - Ellipsoid
0455         \par Projection parameters
0456          - lat_ts: Latitude of true scale (degrees)
0457         \par Example
0458         \image html ex_stere.gif
0459     */
0460     template <typename T, typename Parameters>
0461     struct stere_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
0462     {
0463         template <typename Params>
0464         inline stere_spheroid(Params const& params, Parameters const& par)
0465         {
0466             detail::stere::setup_stere(params, par, this->m_proj_parm);
0467         }
0468     };
0469 
0470     /*!
0471         \brief Universal Polar Stereographic projection
0472         \ingroup projections
0473         \tparam Geographic latlong point type
0474         \tparam Cartesian xy point type
0475         \tparam Parameters parameter type
0476         \par Projection characteristics
0477          - Azimuthal
0478          - Spheroid
0479          - Ellipsoid
0480         \par Projection parameters
0481          - south: Denotes southern hemisphere UTM zone (boolean)
0482         \par Example
0483         \image html ex_ups.gif
0484     */
0485     template <typename T, typename Parameters>
0486     struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
0487     {
0488         template <typename Params>
0489         inline ups_ellipsoid(Params const& params, Parameters & par)
0490         {
0491             detail::stere::setup_ups(params, par, this->m_proj_parm);
0492         }
0493     };
0494 
0495     /*!
0496         \brief Universal Polar Stereographic projection
0497         \ingroup projections
0498         \tparam Geographic latlong point type
0499         \tparam Cartesian xy point type
0500         \tparam Parameters parameter type
0501         \par Projection characteristics
0502          - Azimuthal
0503          - Spheroid
0504          - Ellipsoid
0505         \par Projection parameters
0506          - south: Denotes southern hemisphere UTM zone (boolean)
0507         \par Example
0508         \image html ex_ups.gif
0509     */
0510     template <typename T, typename Parameters>
0511     struct ups_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
0512     {
0513         template <typename Params>
0514         inline ups_spheroid(Params const& params, Parameters & par)
0515         {
0516             detail::stere::setup_ups(params, par, this->m_proj_parm);
0517         }
0518     };
0519 
0520     #ifndef DOXYGEN_NO_DETAIL
0521     namespace detail
0522     {
0523 
0524         // Static projection
0525         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_stere, stere_spheroid, stere_ellipsoid)
0526         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_ups, ups_spheroid, ups_ellipsoid)
0527 
0528         // Factory entry(s)
0529         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(stere_entry, stere_spheroid, stere_ellipsoid)
0530         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(ups_entry, ups_spheroid, ups_ellipsoid)
0531 
0532         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(stere_init)
0533         {
0534             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(stere, stere_entry)
0535             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ups, ups_entry)
0536         }
0537 
0538     } // namespace detail
0539     #endif // doxygen
0540 
0541 } // namespace projections
0542 
0543 }} // namespace boost::geometry
0544 
0545 #endif // BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
0546
[ 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!