[ceph.git] / ceph / src / boost / boost / geometry / formulas / differential_quantities.hpp

// Boost.Geometry

// Copyright (c) 2016-2019 Oracle and/or its affiliates.

// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
#define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP

#include <boost/geometry/core/assert.hpp>

#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>


namespace boost { namespace geometry { namespace formula
{

/*!
\brief The solution of a part of the inverse problem - differential quantities.
\author See
- Charles F.F Karney, Algorithms for geodesics, 2011
https://arxiv.org/pdf/1109.4448.pdf
*/
template <
    typename CT,
    bool EnableReducedLength,
    bool EnableGeodesicScale,
    unsigned int Order = 2,
    bool ApproxF = true
>
class differential_quantities
{
public:
    static inline void apply(CT const& lon1, CT const& lat1,
                             CT const& lon2, CT const& lat2,
                             CT const& azimuth, CT const& reverse_azimuth,
                             CT const& b, CT const& f,
                             CT & reduced_length, CT & geodesic_scale)
    {
        CT const dlon = lon2 - lon1;
        CT const sin_lat1 = sin(lat1);
        CT const cos_lat1 = cos(lat1);
        CT const sin_lat2 = sin(lat2);
        CT const cos_lat2 = cos(lat2);

        apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
              azimuth, reverse_azimuth,
              b, f,
              reduced_length, geodesic_scale);
    }

    static inline void apply(CT const& dlon,
                             CT const& sin_lat1, CT const& cos_lat1,
                             CT const& sin_lat2, CT const& cos_lat2,
                             CT const& azimuth, CT const& reverse_azimuth,
                             CT const& b, CT const& f,
                             CT & reduced_length, CT & geodesic_scale)
    {
        CT const c0 = 0;
        CT const c1 = 1;
        CT const one_minus_f = c1 - f;

        CT sin_bet1 = one_minus_f * sin_lat1;
        CT sin_bet2 = one_minus_f * sin_lat2;
            
        // equator
        if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
        {
            CT const sig_12 = dlon / one_minus_f;
            if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
            {
                BOOST_GEOMETRY_ASSERT((-math::pi<CT>() <= azimuth && azimuth <= math::pi<CT>()));

                int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal
                CT m12 = azi_sign * sin(sig_12) * b;
                reduced_length = m12;
            }
                
            if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
            {
                CT M12 = cos(sig_12);
                geodesic_scale = M12;
            }
        }
        else
        {
            CT const c2 = 2;
            CT const e2 = f * (c2 - f);
            CT const ep2 = e2 / math::sqr(one_minus_f);

            CT const sin_alp1 = sin(azimuth);
            CT const cos_alp1 = cos(azimuth);
            //CT const sin_alp2 = sin(reverse_azimuth);
            CT const cos_alp2 = cos(reverse_azimuth);

            CT cos_bet1 = cos_lat1;
            CT cos_bet2 = cos_lat2;

            normalize(sin_bet1, cos_bet1);
            normalize(sin_bet2, cos_bet2);

            CT sin_sig1 = sin_bet1;
            CT cos_sig1 = cos_alp1 * cos_bet1;
            CT sin_sig2 = sin_bet2;
            CT cos_sig2 = cos_alp2 * cos_bet2;

            normalize(sin_sig1, cos_sig1);
            normalize(sin_sig2, cos_sig2);

            CT const sin_alp0 = sin_alp1 * cos_bet1;
            CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0);

            CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ?
                           J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
                           J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;

            CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
            CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));

            if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
            {
                CT const m12_b = dn2 * (cos_sig1 * sin_sig2)
                               - dn1 * (sin_sig1 * cos_sig2)
                               - cos_sig1 * cos_sig2 * J12;
                CT const m12 = m12_b * b;

                reduced_length = m12;
            }

            if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
            {
                CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2;
                CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2);
                CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1;

                geodesic_scale = M12;
            }
        }
    }

private:
    /*! Approximation of J12, expanded into taylor series in f
        Maxima script:
        ep2: f * (2-f) / ((1-f)^2);
        k2: ca02 * ep2;
        assume(f < 1);
        assume(sig > 0);
        I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
        I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
        J(sig):= I1(sig) - I2(sig);
        S: taylor(J(sig), f, 0, 3);
        S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f );
        S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 );
        S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 );
    */
    static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1,
                           CT const& sin_sig2, CT const& cos_sig2,
                           CT const& cos_alp0_sqr, CT const& f)
    {
        if (Order == 0)
        {
            return 0;
        }

        CT const c2 = 2;

        CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
                                cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2);
        CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
        CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
        CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
        CT const L1 = sig_12 - sin_2sig_12 / c2;

        if (Order == 1)
        {
            return cos_alp0_sqr * f * L1;
        }
        
        CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
        CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
        CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
        
        CT const c8 = 8;
        CT const c12 = 12;
        CT const c16 = 16;
        CT const c24 = 24;

        CT const L2 = -( cos_alp0_sqr * sin_4sig_12
                         + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12
                         + (c12 * cos_alp0_sqr - c24) * sig_12)
                       / c16;

        if (Order == 2)
        {
            return cos_alp0_sqr * f * (L1 + f * L2);
        }

        CT const c4 = 4;
        CT const c9 = 9;
        CT const c48 = 48;
        CT const c60 = 60;
        CT const c64 = 64;
        CT const c96 = 96;
        CT const c128 = 128;
        CT const c144 = 144;

        CT const cos_alp0_quad = math::sqr(cos_alp0_sqr);
        CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
        CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
        CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;

        CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12;
        CT const B = c4 * cos_alp0_quad * sin3_2sig_12;
        CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12;
        CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12;

        CT const L3 = (A + B + C + D) / c64;

        // Order 3 and higher
        return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3));
    }

    /*! Approximation of J12, expanded into taylor series in e'^2
        Maxima script:
        k2: ca02 * ep2;
        assume(sig > 0);
        I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
        I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
        J(sig):= I1(sig) - I2(sig);
        S: taylor(J(sig), ep2, 0, 3);
        S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 );
        S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 );
        S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 );
    */
    static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1,
                                CT const& sin_sig2, CT const& cos_sig2,
                                CT const& cos_alp0_sqr, CT const& ep_sqr)
    {
        if (Order == 0)
        {
            return 0;
        }

        CT const c2 = 2;
        CT const c4 = 4;

        CT const c2a0ep2 = cos_alp0_sqr * ep_sqr;

        CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
                                cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1
        CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
        CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
        CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;

        CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;

        if (Order == 1)
        {
            return c2a0ep2 * L1;
        }

        CT const c8 = 8;
        CT const c64 = 64;
        
        CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
        CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
        CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
        
        CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64;

        if (Order == 2)
        {
            return c2a0ep2 * (L1 + c2a0ep2 * L2);
        }

        CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
        CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
        CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;

        CT const c9 = 9;
        CT const c48 = 48;
        CT const c60 = 60;
        CT const c512 = 512;

        CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512;

        // Order 3 and higher
        return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3));
    }

    static inline void normalize(CT & x, CT & y)
    {
        CT const len = math::sqrt(math::sqr(x) + math::sqr(y));
        x /= len;
        y /= len;
    }
};

}}} // namespace boost::geometry::formula


#endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
Commit	Line	Data
b32b8144 FG	1	// Boost.Geometry
b32b8144 FG	2
92f5a8d4	3	// Copyright (c) 2016-2019 Oracle and/or its affiliates.
b32b8144 FG	4
	5	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	6
	7	// Use, modification and distribution is subject to the Boost Software License,
	8	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	9	// http://www.boost.org/LICENSE_1_0.txt)
	10
	11	#ifndef BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
	12	#define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
	13
92f5a8d4	14	#include <boost/geometry/core/assert.hpp>
b32b8144 FG	15
	16	#include <boost/geometry/util/condition.hpp>
	17	#include <boost/geometry/util/math.hpp>
	18
	19
	20	namespace boost { namespace geometry { namespace formula
	21	{
	22
	23	/*!
	24	\brief The solution of a part of the inverse problem - differential quantities.
	25	\author See
	26	- Charles F.F Karney, Algorithms for geodesics, 2011
	27	https://arxiv.org/pdf/1109.4448.pdf
	28	*/
	29	template <
	30	typename CT,
	31	bool EnableReducedLength,
	32	bool EnableGeodesicScale,
	33	unsigned int Order = 2,
	34	bool ApproxF = true
	35	>
	36	class differential_quantities
	37	{
	38	public:
	39	static inline void apply(CT const& lon1, CT const& lat1,
	40	CT const& lon2, CT const& lat2,
	41	CT const& azimuth, CT const& reverse_azimuth,
	42	CT const& b, CT const& f,
	43	CT & reduced_length, CT & geodesic_scale)
	44	{
	45	CT const dlon = lon2 - lon1;
	46	CT const sin_lat1 = sin(lat1);
	47	CT const cos_lat1 = cos(lat1);
	48	CT const sin_lat2 = sin(lat2);
	49	CT const cos_lat2 = cos(lat2);
	50
	51	apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
	52	azimuth, reverse_azimuth,
	53	b, f,
	54	reduced_length, geodesic_scale);
	55	}
	56
	57	static inline void apply(CT const& dlon,
	58	CT const& sin_lat1, CT const& cos_lat1,
	59	CT const& sin_lat2, CT const& cos_lat2,
	60	CT const& azimuth, CT const& reverse_azimuth,
	61	CT const& b, CT const& f,
	62	CT & reduced_length, CT & geodesic_scale)
	63	{
	64	CT const c0 = 0;
	65	CT const c1 = 1;
	66	CT const one_minus_f = c1 - f;
	67
	68	CT sin_bet1 = one_minus_f * sin_lat1;
	69	CT sin_bet2 = one_minus_f * sin_lat2;
	70
	71	// equator
	72	if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
	73	{
92f5a8d4	74	CT const sig_12 = dlon / one_minus_f;
b32b8144 FG	75	if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
b32b8144 FG	76	{
92f5a8d4 TL	77	BOOST_GEOMETRY_ASSERT((-math::pi<CT>() <= azimuth && azimuth <= math::pi<CT>()));
	78
	79	int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal
	80	CT m12 = azi_sign * sin(sig_12) * b;
b32b8144 FG	81	reduced_length = m12;
	82	}
	83
	84	if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
	85	{
	86	CT M12 = cos(sig_12);
	87	geodesic_scale = M12;
	88	}
	89	}
	90	else
	91	{
	92	CT const c2 = 2;
	93	CT const e2 = f * (c2 - f);
	94	CT const ep2 = e2 / math::sqr(one_minus_f);
	95
	96	CT const sin_alp1 = sin(azimuth);
	97	CT const cos_alp1 = cos(azimuth);
	98	//CT const sin_alp2 = sin(reverse_azimuth);
	99	CT const cos_alp2 = cos(reverse_azimuth);
	100
	101	CT cos_bet1 = cos_lat1;
	102	CT cos_bet2 = cos_lat2;
	103
	104	normalize(sin_bet1, cos_bet1);
	105	normalize(sin_bet2, cos_bet2);
	106
	107	CT sin_sig1 = sin_bet1;
	108	CT cos_sig1 = cos_alp1 * cos_bet1;
	109	CT sin_sig2 = sin_bet2;
	110	CT cos_sig2 = cos_alp2 * cos_bet2;
	111
	112	normalize(sin_sig1, cos_sig1);
	113	normalize(sin_sig2, cos_sig2);
	114
	115	CT const sin_alp0 = sin_alp1 * cos_bet1;
	116	CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0);
	117
	118	CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ?
	119	J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
	120	J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;
	121
	122	CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
	123	CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
	124
	125	if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
	126	{
	127	CT const m12_b = dn2 * (cos_sig1 * sin_sig2)
	128	- dn1 * (sin_sig1 * cos_sig2)
	129	- cos_sig1 * cos_sig2 * J12;
	130	CT const m12 = m12_b * b;
	131
	132	reduced_length = m12;
	133	}
	134
	135	if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
	136	{
	137	CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2;
	138	CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2);
	139	CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1;
	140
	141	geodesic_scale = M12;
	142	}
	143	}
	144	}
145
146	private:
147	/*! Approximation of J12, expanded into taylor series in f
148	Maxima script:
149	ep2: f * (2-f) / ((1-f)^2);
150	k2: ca02 * ep2;
151	assume(f < 1);
152	assume(sig > 0);
153	I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
154	I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
155	J(sig):= I1(sig) - I2(sig);
156	S: taylor(J(sig), f, 0, 3);
157	S1: factor( 2integrate(sin(s)^2,s,0,sig)ca02*f );
158	S2: factor( ((integrate(-6ca02^2sin(s)^4+6ca02sin(s)^2,s,0,sig)+integrate(-2ca02^2sin(s)^4+6ca02sin(s)^2,s,0,sig))*f^2)/4 );
159	S3: factor( ((integrate(30ca02^3sin(s)^6-54ca02^2sin(s)^4+24ca02sin(s)^2,s,0,sig)+integrate(6ca02^3sin(s)^6-18ca02^2sin(s)^4+24ca02sin(s)^2,s,0,sig))*f^3)/12 );
160	*/
161	static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1,
162	CT const& sin_sig2, CT const& cos_sig2,
163	CT const& cos_alp0_sqr, CT const& f)
164	{
165	if (Order == 0)
166	{
167	return 0;
168	}
169
170	CT const c2 = 2;
171
172	CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
173	cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2);
174	CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
175	CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
176	CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
177	CT const L1 = sig_12 - sin_2sig_12 / c2;
178
179	if (Order == 1)
180	{
181	return cos_alp0_sqr * f * L1;
182	}
183
184	CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
185	CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
186	CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
187
188	CT const c8 = 8;
189	CT const c12 = 12;
190	CT const c16 = 16;
191	CT const c24 = 24;
192
193	CT const L2 = -( cos_alp0_sqr * sin_4sig_12
194	+ (-c8 * cos_alp0_sqr + c12) * sin_2sig_12
195	+ (c12 * cos_alp0_sqr - c24) * sig_12)
196	/ c16;
197
198	if (Order == 2)
199	{
200	return cos_alp0_sqr * f * (L1 + f * L2);
201	}
202
203	CT const c4 = 4;
204	CT const c9 = 9;
205	CT const c48 = 48;
206	CT const c60 = 60;
207	CT const c64 = 64;
208	CT const c96 = 96;
209	CT const c128 = 128;
210	CT const c144 = 144;
211
212	CT const cos_alp0_quad = math::sqr(cos_alp0_sqr);
213	CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
214	CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
215	CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
216
217	CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12;
218	CT const B = c4 * cos_alp0_quad * sin3_2sig_12;
219	CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12;
220	CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12;
221
222	CT const L3 = (A + B + C + D) / c64;
223
224	// Order 3 and higher
225	return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3));
226	}
227
228	/*! Approximation of J12, expanded into taylor series in e'^2
229	Maxima script:
230	k2: ca02 * ep2;
231	assume(sig > 0);
232	I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
233	I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
234	J(sig):= I1(sig) - I2(sig);
235	S: taylor(J(sig), ep2, 0, 3);
236	S1: factor( integrate(sin(s)^2,s,0,sig)ca02ep2 );
237	S2: factor( (integrate(sin(s)^4,s,0,sig)ca02^2ep2^2)/2 );
238	S3: factor( (3integrate(sin(s)^6,s,0,sig)ca02^3*ep2^3)/8 );
239	*/
240	static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1,
241	CT const& sin_sig2, CT const& cos_sig2,
242	CT const& cos_alp0_sqr, CT const& ep_sqr)
243	{
244	if (Order == 0)
245	{
246	return 0;
247	}
248
249	CT const c2 = 2;
250	CT const c4 = 4;
251
252	CT const c2a0ep2 = cos_alp0_sqr * ep_sqr;
253
254	CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
255	cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1
256	CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
257	CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
258	CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
259
260	CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;
261
262	if (Order == 1)
263	{
264	return c2a0ep2 * L1;
265	}
266
267	CT const c8 = 8;
268	CT const c64 = 64;
269
270	CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
271	CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
272	CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
273
274	CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64;
275
276	if (Order == 2)
277	{
278	return c2a0ep2 * (L1 + c2a0ep2 * L2);
279	}
280
281	CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
282	CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
283	CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
284
285	CT const c9 = 9;
286	CT const c48 = 48;
287	CT const c60 = 60;
288	CT const c512 = 512;
289
290	CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512;
291
292	// Order 3 and higher
293	return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3));
294	}
295
296	static inline void normalize(CT & x, CT & y)
297	{
298	CT const len = math::sqrt(math::sqr(x) + math::sqr(y));
299	x /= len;
300	y /= len;
301	}
302	};
303
304	}}} // namespace boost::geometry::formula
305
306
307	#endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP