Ticket #12185: transverse-mercator.patch

File transverse-mercator.patch, 22.5 KB (added by bastiK, 9 years ago)

  • src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection.proj;
     3
     4import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
     5
     6/**
     7 * Abstract base class providing utilities for implementations of the Proj
     8 * interface.
     9 *
     10 * This class has been derived from the implementation of the Geotools project;
     11 * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
     12 * at the time of migration.
     13 * <p>
     14 *
     15 * @author André Gosselin
     16 * @author Martin Desruisseaux (PMO, IRD)
     17 * @author Rueben Schulz
     18*/
     19public abstract class AbstractProj implements Proj {
     20
     21    /**
     22     * Maximum number of iterations for iterative computations.
     23     */
     24    private static final int MAXIMUM_ITERATIONS = 15;
     25
     26    /**
     27     * Relative iteration precision used in the <code>mlfn<code> method
     28     */
     29    private static final double MLFN_TOL = 1E-11;
     30
     31    /**
     32     * Constants used to calculate {@link #en0}, {@link #en1},
     33     * {@link #en2}, {@link #en3}, {@link #en4}.
     34     */
     35    private static final double C00= 1.0,
     36                                C02= 0.25,
     37                                C04= 0.046875,
     38                                C06= 0.01953125,
     39                                C08= 0.01068115234375,
     40                                C22= 0.75,
     41                                C44= 0.46875,
     42                                C46= 0.01302083333333333333,
     43                                C48= 0.00712076822916666666,
     44                                C66= 0.36458333333333333333,
     45                                C68= 0.00569661458333333333,
     46                                C88= 0.3076171875;
     47
     48    /**
     49     * Constant needed for the <code>mlfn<code> method.
     50     * Setup at construction time.
     51     */
     52    protected double en0,en1,en2,en3,en4;
     53
     54    /**
     55     * The square of excentricity: e² = (a²-b²)/a² where
     56     * <var>e</var> is the {@linkplain #excentricity excentricity},
     57     * <var>a</var> is the {@linkplain #semiMajor semi major} axis length and
     58     * <var>b</var> is the {@linkplain #semiMinor semi minor} axis length.
     59     */
     60    protected double e2;
     61
     62    @Override
     63    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
     64        e2 = params.ellps.e2;
     65        //  Compute constants for the mlfn
     66        double t;
     67        en0 = C00 - e2  *  (C02 + e2  *
     68             (C04 + e2  *  (C06 + e2  * C08)));
     69        en1 =       e2  *  (C22 - e2  *
     70             (C04 + e2  *  (C06 + e2  * C08)));
     71        en2 =  (t = e2  *         e2) *
     72             (C44 - e2  *  (C46 + e2  * C48));
     73        en3 = (t *= e2) *  (C66 - e2  * C68);
     74        en4 =   t * e2  *  C88;
     75    }
     76
     77    /**
     78     * Calculates the meridian distance. This is the distance along the central
     79     * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters
     80     * when used in conjuction with typical major axis values.
     81     *
     82     * @param phi latitude to calculate meridian distance for.
     83     * @param sphi sin(phi).
     84     * @param cphi cos(phi).
     85     * @return meridian distance for the given latitude.
     86     */
     87    protected final double mlfn(final double phi, double sphi, double cphi) {
     88        cphi *= sphi;
     89        sphi *= sphi;
     90        return en0 * phi - cphi *
     91              (en1 + sphi *
     92              (en2 + sphi *
     93              (en3 + sphi *
     94              (en4))));
     95    }
     96
     97    /**
     98     * Calculates the latitude ({@code phi}) from a meridian distance.
     99     * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
     100     *
     101     * @param arg meridian distance to calulate latitude for.
     102     * @return the latitude of the meridian distance.
     103     * @throws RuntimeException if the itteration does not converge.
     104     */
     105    protected final double inv_mlfn(double arg) {
     106        double s, t, phi, k = 1.0/(1.0 - e2);
     107        int i;
     108        phi = arg;
     109        for (i=MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
     110            if (--i < 0) {
     111                throw new RuntimeException("Too many iterations");
     112            }
     113            s = Math.sin(phi);
     114            t = 1.0 - e2 * s * s;
     115            t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
     116            phi -= t;
     117            if (Math.abs(t) < MLFN_TOL) {
     118                return phi;
     119            }
     120        }
     121    }
     122
     123}

  • src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java 

    Property changes on: src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
     
    11// License: GPL. For details, see LICENSE file.
    22package org.openstreetmap.josm.data.projection.proj;
    33
    4 import static java.lang.Math.cos;
    5 import static java.lang.Math.pow;
    6 import static java.lang.Math.sin;
    7 import static java.lang.Math.sqrt;
    8 import static java.lang.Math.tan;
    94import static org.openstreetmap.josm.tools.I18n.tr;
    105
    116import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
    127
    138/**
    14  * Transverse Mercator projection.
     9 * Transverse Mercator Projection (EPSG code 9807). This
     10 * is a cylindrical projection, in which the cylinder has been rotated 90°.
     11 * Instead of being tangent to the equator (or to an other standard latitude),
     12 * it is tangent to a central meridian. Deformation are more important as we
     13 * are going futher from the central meridian. The Transverse Mercator
     14 * projection is appropriate for region wich have a greater extent north-south
     15 * than east-west.
     16 * <p>
    1517 *
    16  * @author Dirk Stöcker
    17  * code based on JavaScript from Chuck Taylor
     18 * The elliptical equations used here are series approximations, and their accuracy
     19 * decreases as points move farther from the central meridian of the projection.
     20 * The forward equations here are accurate to a less than a mm &plusmn;10 degrees from
     21 * the central meridian, a few mm &plusmn;15 degrees from the
     22 * central meridian and a few cm &plusmn;20 degrees from the central meridian.
     23 * The spherical equations are not approximations and should always give the
     24 * correct values.
     25 * <p>
    1826 *
     27 * There are a number of versions of the transverse mercator projection
     28 * including the Universal (UTM) and Modified (MTM) Transverses Mercator
     29 * projections. In these cases the earth is divided into zones. For the UTM
     30 * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
     31 * 180 degrees longitude, and between lats 84 degrees North and 80
     32 * degrees South. The central meridian is taken as the center of the zone
     33 * and the latitude of origin is the equator. A scale factor of 0.9996 and
     34 * false easting of 500000m is used for all zones and a false northing of 10000000m
     35 * is used for zones in the southern hemisphere.
     36 * <p>
     37 *
     38 * NOTE: formulas used below are not those of Snyder, but rather those
     39 *       from the {@code proj4} package of the USGS survey, which
     40 *       have been reproduced verbatim. USGS work is acknowledged here.
     41 * <p>
     42 *
     43 * This class has been derived from the implementation of the Geotools project;
     44 * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator
     45 * at the time of migration.
     46 * <p>
     47 *
     48 * <b>References:</b>
     49 * <ul>
     50 *   <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
     51 *        Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
     52 *   <li> John P. Snyder (Map Projections - A Working Manual,
     53 *        U.S. Geological Survey Professional Paper 1395, 1987).</li>
     54 *   <li> "Coordinate Conversions and Transformations including Formulas",
     55 *        EPSG Guidence Note Number 7, Version 19.</li>
     56 * </ul>
     57 *
     58 * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
     59 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
     60 *
     61 * @author André Gosselin
     62 * @author Martin Desruisseaux (PMO, IRD)
     63 * @author Rueben Schulz
    1964 */
    20 public class TransverseMercator implements Proj {
     65public class TransverseMercator extends AbstractProj {
    2166
    22     protected double a, b;
     67    /**
     68     * Contants used for the forward and inverse transform for the eliptical
     69     * case of the Transverse Mercator.
     70     */
     71    private static final double FC1= 1.00000000000000000000000,  // 1/1
     72                                FC2= 0.50000000000000000000000,  // 1/2
     73                                FC3= 0.16666666666666666666666,  // 1/6
     74                                FC4= 0.08333333333333333333333,  // 1/12
     75                                FC5= 0.05000000000000000000000,  // 1/20
     76                                FC6= 0.03333333333333333333333,  // 1/30
     77                                FC7= 0.02380952380952380952380,  // 1/42
     78                                FC8= 0.01785714285714285714285;  // 1/56
    2379
     80    /**
     81     * Maximum difference allowed when comparing real numbers.
     82     */
     83    private static final double EPSILON = 1E-6;
     84
     85    /**
     86     * A derived quantity of excentricity, computed by <code>e'² = (a²-b²)/b² = es/(1-es)</code>
     87     * where <var>a</var> is the semi-major axis length and <var>b</bar> is the semi-minor axis
     88     * length.
     89     */
     90    private double eb2;
     91
     92    /**
     93     * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
     94     * This is called '<var>phi0</var>' in Snyder.
     95     * <p>
     96     * <strong>Consider this field as final</strong>. It is not final only
     97     * because some classes need to modify it at construction time.
     98     */
     99    protected double latitudeOfOrigin;
     100
     101    /**
     102     * Meridian distance at the {@code latitudeOfOrigin}.
     103     * Used for calculations for the ellipsoid.
     104     */
     105    private double ml0;
     106
    24107    @Override
    25108    public String getName() {
    26109        return tr("Transverse Mercator");
    … …  
    33116
    34117    @Override
    35118    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
    36         this.a = params.ellps.a;
    37         this.b = params.ellps.b;
     119        super.initialize(params);
     120        eb2 = params.ellps.eb2;
     121        latitudeOfOrigin = params.lat0 == null ? 0 : Math.toRadians(params.lat0);
     122        ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin));
    38123    }
    39124
    40     /**
    41      * Converts a latitude/longitude pair to x and y coordinates in the
    42      * Transverse Mercator projection.  Note that Transverse Mercator is not
    43      * the same as UTM; a scale factor is required to convert between them.
    44      *
    45      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    46      * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    47      *
    48      * @param phi Latitude of the point, in radians
    49      * @param lambda Longitude of the point, in radians
    50      * @return A 2-element array containing the x and y coordinates
    51      *         of the computed point
    52      */
    53125    @Override
    54     public double[] project(double phi, double lambda) {
     126    public double[] project(double y, double x) {
     127        double sinphi = Math.sin(y);
     128        double cosphi = Math.cos(y);
    55129
    56         /* Precalculate ep2 */
    57         double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0);
     130        double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0;
     131        t *= t;
     132        double al = cosphi*x;
     133        double als = al*al;
     134        al /= Math.sqrt(1.0 - e2 * sinphi*sinphi);
     135        double n = eb2 * cosphi*cosphi;
    58136
    59         /* Precalculate nu2 */
    60         double nu2 = ep2 * pow(cos(phi), 2.0);
     137        /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
     138        y = (mlfn(y, sinphi, cosphi) - ml0 +
     139            sinphi * al * x *
     140            FC2 * ( 1.0 +
     141            FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) +
     142            FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) +
     143            FC8 * als * (1385.0 + t * ( t*(543.0 - t) - 3111.0))))));
    61144
    62         /* Precalculate N / a */
    63         double N_a = a / (b * sqrt(1 + nu2));
     145        x = al*(FC1 + FC3 * als*(1.0 - t + n +
     146            FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) +
     147            FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0 )))));
    64148
    65         /* Precalculate t */
    66         double t = tan(phi);
    67         double t2 = t * t;
    68 
    69         /* Precalculate l */
    70         double l = lambda;
    71 
    72         /* Precalculate coefficients for l**n in the equations below
    73            so a normal human being can read the expressions for easting
    74            and northing
    75            -- l**1 and l**2 have coefficients of 1.0 */
    76         double l3coef = 1.0 - t2 + nu2;
    77 
    78         double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
    79 
    80         double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
    81         - 58.0 * t2 * nu2;
    82 
    83         double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
    84         - 330.0 * t2 * nu2;
    85 
    86         double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
    87 
    88         double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
    89 
    90         return new double[] {
    91                 /* Calculate easting (x) */
    92                 N_a * cos(phi) * l
    93                 + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
    94                 + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
    95                 + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)),
    96                 /* Calculate northing (y) */
    97                 ArcLengthOfMeridian(phi) / a
    98                 + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0))
    99                 + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
    100                 + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
    101                 + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) };
     149        return new double[] { x, y };
    102150    }
    103151
    104     /**
    105      * Converts x and y coordinates in the Transverse Mercator projection to
    106      * a latitude/longitude pair.  Note that Transverse Mercator is not
    107      * the same as UTM; a scale factor is required to convert between them.
    108      *
    109      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    110      *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    111      *
    112      * Remarks:
    113      *   The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
    114      *   N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
    115      *   to the footpoint latitude phif.
    116      *
    117      *   x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
    118      *   to optimize computations.
    119      *
    120      * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid
    121      * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid
    122      * @return A 2-element containing the latitude and longitude
    123      *               in radians
    124      */
    125152    @Override
    126153    public double[] invproject(double x, double y) {
    127         /* Get the value of phif, the footpoint latitude. */
    128         double phif = footpointLatitude(y);
     154        double phi = inv_mlfn(ml0 + y);
    129155
    130         /* Precalculate ep2 */
    131         double ep2 = (a*a - b*b)
    132         / (b*b);
     156        if (Math.abs(phi) >= Math.PI/2) {
     157            y = y<0.0 ? -(Math.PI/2) : (Math.PI/2);
     158            x = 0.0;
     159        } else {
     160            double sinphi = Math.sin(phi);
     161            double cosphi = Math.cos(phi);
     162            double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0;
     163            double n = eb2 * cosphi*cosphi;
     164            double con = 1.0 - e2 * sinphi*sinphi;
     165            double d = x * Math.sqrt(con);
     166            con *= t;
     167            t *= t;
     168            double ds = d*d;
    133169
    134         /* Precalculate cos(phif) */
    135         double cf = cos(phif);
     170            y = phi - (con*ds / (1.0 - e2)) *
     171                FC2 * (1.0 - ds *
     172                FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds *
     173                FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds *
     174                FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t))))));
    136175
    137         /* Precalculate nuf2 */
    138         double nuf2 = ep2 * pow(cf, 2.0);
    139 
    140         /* Precalculate Nf / a and initialize Nfpow */
    141         double Nf_a = a / (b * sqrt(1 + nuf2));
    142         double Nfpow = Nf_a;
    143 
    144         /* Precalculate tf */
    145         double tf = tan(phif);
    146         double tf2 = tf * tf;
    147         double tf4 = tf2 * tf2;
    148 
    149         /* Precalculate fractional coefficients for x**n in the equations
    150            below to simplify the expressions for latitude and longitude. */
    151         double x1frac = 1.0 / (Nfpow * cf);
    152 
    153         Nfpow *= Nf_a;   /* now equals Nf**2) */
    154         double x2frac = tf / (2.0 * Nfpow);
    155 
    156         Nfpow *= Nf_a;   /* now equals Nf**3) */
    157         double x3frac = 1.0 / (6.0 * Nfpow * cf);
    158 
    159         Nfpow *= Nf_a;   /* now equals Nf**4) */
    160         double x4frac = tf / (24.0 * Nfpow);
    161 
    162         Nfpow *= Nf_a;   /* now equals Nf**5) */
    163         double x5frac = 1.0 / (120.0 * Nfpow * cf);
    164 
    165         Nfpow *= Nf_a;   /* now equals Nf**6) */
    166         double x6frac = tf / (720.0 * Nfpow);
    167 
    168         Nfpow *= Nf_a;   /* now equals Nf**7) */
    169         double x7frac = 1.0 / (5040.0 * Nfpow * cf);
    170 
    171         Nfpow *= Nf_a;   /* now equals Nf**8) */
    172         double x8frac = tf / (40320.0 * Nfpow);
    173 
    174         /* Precalculate polynomial coefficients for x**n.
    175            -- x**1 does not have a polynomial coefficient. */
    176         double x2poly = -1.0 - nuf2;
    177         double x3poly = -1.0 - 2 * tf2 - nuf2;
    178         double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
    179         double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
    180         double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
    181         double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
    182         double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
    183 
    184         return new double[] {
    185                         /* Calculate latitude */
    186                         phif + x2frac * x2poly * (x * x)
    187                         + x4frac * x4poly * pow(x, 4.0)
    188                         + x6frac * x6poly * pow(x, 6.0)
    189                         + x8frac * x8poly * pow(x, 8.0),
    190                         /* Calculate longitude */
    191                         x1frac * x
    192                         + x3frac * x3poly * pow(x, 3.0)
    193                         + x5frac * x5poly * pow(x, 5.0)
    194                         + x7frac * x7poly * pow(x, 7.0) };
     176            x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n -
     177                ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n -
     178                ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi;
     179        }
     180        return new double[] { y, x };
    195181    }
    196 
    197     /**
    198      * ArcLengthOfMeridian
    199      *
    200      * Computes the ellipsoidal distance from the equator to a point at a
    201      * given latitude.
    202      *
    203      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    204      * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    205      *
    206      * @param phi Latitude of the point, in radians
    207      * @return The ellipsoidal distance of the point from the equator
    208      *         (in meters, divided by the semi major axis of the ellipsoid)
    209      */
    210     private double ArcLengthOfMeridian(double phi) {
    211         /* Precalculate n */
    212         double n = (a - b) / (a + b);
    213 
    214         /* Precalculate alpha */
    215         double alpha = ((a + b) / 2.0)
    216             * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
    217 
    218         /* Precalculate beta */
    219         double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
    220             + (-3.0 * pow(n, 5.0) / 32.0);
    221 
    222         /* Precalculate gamma */
    223         double gamma = (15.0 * pow(n, 2.0) / 16.0)
    224             + (-15.0 * pow(n, 4.0) / 32.0);
    225 
    226         /* Precalculate delta */
    227         double delta = (-35.0 * pow(n, 3.0) / 48.0)
    228             + (105.0 * pow(n, 5.0) / 256.0);
    229 
    230         /* Precalculate epsilon */
    231         double epsilon = 315.0 * pow(n, 4.0) / 512.0;
    232 
    233         /* Now calculate the sum of the series and return */
    234         return alpha
    235             * (phi + (beta * sin(2.0 * phi))
    236                     + (gamma * sin(4.0 * phi))
    237                     + (delta * sin(6.0 * phi))
    238                     + (epsilon * sin(8.0 * phi)));
    239     }
    240 
    241     /**
    242      * FootpointLatitude
    243      *
    244      * Computes the footpoint latitude for use in converting transverse
    245      * Mercator coordinates to ellipsoidal coordinates.
    246      *
    247      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    248      *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    249      *
    250      * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid
    251      * @return The footpoint latitude, in radians
    252      */
    253     private double footpointLatitude(double y) {
    254         /* Precalculate n (Eq. 10.18) */
    255         double n = (a - b) / (a + b);
    256 
    257         /* Precalculate alpha_ (Eq. 10.22) */
    258         /* (Same as alpha in Eq. 10.17) */
    259         double alpha_ = ((a + b) / 2.0)
    260             * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
    261 
    262         /* Precalculate y_ (Eq. 10.23) */
    263         double y_ = y / alpha_ * a;
    264 
    265         /* Precalculate beta_ (Eq. 10.22) */
    266         double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
    267             + (269.0 * pow(n, 5.0) / 512.0);
    268 
    269         /* Precalculate gamma_ (Eq. 10.22) */
    270         double gamma_ = (21.0 * pow(n, 2.0) / 16.0)
    271             + (-55.0 * pow(n, 4.0) / 32.0);
    272 
    273         /* Precalculate delta_ (Eq. 10.22) */
    274         double delta_ = (151.0 * pow(n, 3.0) / 96.0)
    275             + (-417.0 * pow(n, 5.0) / 128.0);
    276 
    277         /* Precalculate epsilon_ (Eq. 10.22) */
    278         double epsilon_ = 1097.0 * pow(n, 4.0) / 512.0;
    279 
    280         /* Now calculate the sum of the series (Eq. 10.21) */
    281         return y_ + (beta_ * sin(2.0 * y_))
    282             + (gamma_ * sin(4.0 * y_))
    283             + (delta_ * sin(6.0 * y_))
    284             + (epsilon_ * sin(8.0 * y_));
    285     }
    286182}