    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_testing/src/org/deegree/crs/projections/ProjectionUtils.java $
    002    /*----------------------------------------------------------------------------
    003     This file is part of deegree, http://deegree.org/
    004     Copyright (C) 2001-2009 by:
    005       Department of Geography, University of Bonn
    006     and
    007       lat/lon GmbH
    008
    009     This library is free software; you can redistribute it and/or modify it under
    010     the terms of the GNU Lesser General Public License as published by the Free
    011     Software Foundation; either version 2.1 of the License, or (at your option)
    012     any later version.
    013     This library is distributed in the hope that it will be useful, but WITHOUT
    014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
    015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
    016     details.
    017     You should have received a copy of the GNU Lesser General Public License
    018     along with this library; if not, write to the Free Software Foundation, Inc.,
    019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
    020
    021     Contact information:
    022
    023     lat/lon GmbH
    024     Aennchenstr. 19, 53177 Bonn
    025     Germany
    026     http://lat-lon.de/
    027
    028     Department of Geography, University of Bonn
    029     Prof. Dr. Klaus Greve
    030     Postfach 1147, 53001 Bonn
    031     Germany
    032     http://www.geographie.uni-bonn.de/deegree/
    033
    034     e-mail: info@deegree.org
    035    ----------------------------------------------------------------------------*/
    036
    037    package org.deegree.crs.projections;
    038
    039    import java.awt.geom.Rectangle2D;
    040
    041    /**
    042     * The <code>Utils</code> class combines some helpful constants and forms.
    043     *
    044     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
    045     *
    046     * @author last edited by: $Author: mschneider $
    047     *
    048     * @version $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18. Jun 2009) $
    049     *
    050     */
    051
    052    public class ProjectionUtils {
    053
    054        // Very handy for setting maximum number of iterations in a for loop.
    055        private final static int MAX_ITER = 10;
    056
    057        /**
    058         * A small epsilon value
    059         */
    060        public final static double EPS10 = 1e-10;
    061
    062        /**
    063         * An even smaller epsilon value
    064         */
    065        public final static double EPS11 = 1e-11;
    066
    067        /**
    068         * Cotaining the value 0.5*pi
    069         */
    070        public final static double HALFPI = Math.PI * 0.5;
    071
    072        /**
    073         * Containing the value 0.25*pi
    074         */
    075        public final static double QUARTERPI = Math.PI * 0.25;
    076
    077        /**
    078         * Containing the value 2*pi
    079         */
    080        public final static double TWOPI = Math.PI * 2.0;
    081
    082        /**
    083         * Radians to Degrees (180.0/Math.PI)
    084         */
    085        public final static double RTD = 180.0 / Math.PI;
    086
    087        /**
    088         * Degrees to Radians (Math.PI/180.0)
    089         */
    090        public final static double DTR = Math.PI / 180.0;
    091
    092        /**
    093         * The max and min of the projected word map in radians (-Math.PI, -HALFPI, TWOPI, Math.PI)
    094         */
    095        public final static Rectangle2D WORLD_BOUNDS_RAD = new Rectangle2D.Double( -Math.PI, -HALFPI, TWOPI, Math.PI );
    096
    097        /**
    098         * The max and min of the projected word map in degrees (-180, -90, 360, 180)
    099         */
    100        public final static Rectangle2D WORLD_BOUNDS = new Rectangle2D.Double( -180, -90, 360, 180 );
    101
    102        /**
    103         * From the proj4 library, to determine small q which is needed to calculate the authalic (equal-areaed) latitude
    104         * beta, on a sphere having the same surface area as the ellipsoid, relative to the ellipsoid. Snyder (3 -12).
    105         *
    106         * @param sinphi
    107         *            the sine of the angle between the positive z-axis and the line formed between the origin and P.
    108         * @param e
    109         *            the eccentricity
    110         * @return the q value from Snyder (3-12)
    111         * @deprecated use {@link ProjectionUtils#calcQForAuthalicLatitude(double, double)} instead.
    112         */
    113        @Deprecated
    114        public static double qsfn( double sinphi, double e ) {
    115            return calcQForAuthalicLatitude( sinphi, e );
    116        }
    117
    118        /**
    119         * From the proj4 library, to determine small q which is needed to calculate the authalic (equal-areaed) latitude
    120         * beta, on a sphere having the same surface area as the ellipsoid, relative to the ellipsoid. Snyder (3 -12).
    121         *
    122         * @param sinphi
    123         *            the sine of the angle between the positive z-axis and the line formed between the origin and P.
    124         * @param eccentricity
    125         *            the eccentricity of the ellipsoid to map the sphere to.
    126         * @return the q value from Snyder (3-12)
    127         */
    128        public static double calcQForAuthalicLatitude( double sinphi, double eccentricity ) {
    129            if ( eccentricity >= EPS10 ) {
    130                double eAndSinphi = eccentricity * sinphi;
    131                double es = eccentricity * eccentricity;
    132                return ( 1 - es )
    133                       * ( sinphi / ( 1. - eAndSinphi * eAndSinphi ) - ( .5 / eccentricity )
    134                                                                       * Math.log( ( 1. - eAndSinphi ) / ( 1. + eAndSinphi ) ) );
    135            }
    136            // we have a sphere.
    137            return ( sinphi + sinphi );
    138        }
    139
    140        /**
    141         * Pre-calculated values for Snyder's formula (3-5)
    142         */
    143        private final static double T00 = 0.5;
    144
    145        private final static double T01 = .20833333333333333333;/* 5/24. */
    146
    147        private final static double T02 = .08333333333333333333;/* 1/12. */
    148
    149        private final static double T03 = .03611111111111111111;/* 13/360 */
    150
    151        private final static double T10 = .14583333333333333333;/* 7/48 */
    152
    153        private final static double T11 = .12083333333333333333;/* 29/240 */
    154
    155        private final static double T12 = .07039930555555555555;/* 811/11520 */
    156
    157        private final static double T20 = .05833333333333333333;/* 7/120 */
    158
    159        private final static double T21 = .07232142857142857142;/* 81/1120 */
    160
    161        private final static double T30 = .02653149801587301587;/* 4279/161280 */
    162
    163        /**
    164         * Pre-Calculates the values (used for the adams? series) which will be used to calculate the phi value of an
    165         * inverse projection. Snyder (3-5).
    166         *
    167         * @param eccentricitySquared
    168         *            the squared eccentricity from the ellipsoid to calculate the theta for.
    169         * @return the precalculated values.
    170         */
    171        public static double[] preCalcedThetaSeries( double eccentricitySquared ) {
    172            double[] precalculatedSerie = new double[4];
    173            precalculatedSerie[0] = eccentricitySquared * T00;
    174
    175            // eccentricity^4
    176            double tmp = eccentricitySquared * eccentricitySquared;
    177            precalculatedSerie[0] += tmp * T01;
    178            precalculatedSerie[1] = tmp * T10;
    179
    180            // eccentricity^6
    181            tmp *= eccentricitySquared;
    182            precalculatedSerie[0] += tmp * T02;
    183            precalculatedSerie[1] += tmp * T11;
    184            precalculatedSerie[2] = tmp * T20;
    185
    186            // eccentricity^8
    187            tmp *= eccentricitySquared;
    188            precalculatedSerie[0] += tmp * T03;
    189            precalculatedSerie[1] += tmp * T12;
    190            precalculatedSerie[2] += tmp * T21;
    191            precalculatedSerie[3] = tmp * T30;
    192
    193            return precalculatedSerie;
    194        }
    195
    196        /**
    197         * Gets Phi from the given conformal latitude chi and the precalculated values (gotten from
    198         * {@link ProjectionUtils#preCalcedThetaSeries(double)} ) of the adams? serie. From Snyder (3-5).
    199         *
    200         * @param chi
    201         *            the conformal latitude
    202         * @param APA
    203         *            the precalculated values from the serie gotten from
    204         *            {@link ProjectionUtils#preCalcedThetaSeries(double)}.
    205         * @return the Phi as a polarcoordinate on the ellipsoid or chi if the length of APA != 4.
    206         */
    207        public static double calcPhiFromConformalLatitude( double chi, double[] APA ) {
    208            if ( APA.length != 4 ) {
    209                return chi;
    210            }
    211            double tmp = chi + chi;
    212            return ( chi + APA[0] * Math.sin( tmp ) + APA[1] * Math.sin( tmp + tmp ) + APA[2] * Math.sin( tmp + tmp + tmp ) + APA[3]
    213                                                                                                                              * Math.sin( tmp
    214                                                                                                                                          + tmp
    215                                                                                                                                          + tmp
    216                                                                                                                                          + tmp ) );
    217        }
    218
    219        /**
    220         * P[0][0-2] = 1/3, 31/180, 517/5040, P[1][0-2] = 23/360, 251/3780 P[2][0] = 761/45360
    221         */
    222        private final static double P00 = .33333333333333333333; /* 1/3 */
    223
    224        private final static double P01 = .17222222222222222222; /* 31 / 180 */
    225
    226        private final static double P02 = .10257936507936507936; /* 517 / 5040 */
    227
    228        private final static double P10 = .06388888888888888888; /* 23/360 */
    229
    230        private final static double P11 = .06640211640211640211; /* 251/3780 */
    231
    232        private final static double P20 = .01641501294219154443; /* 761/45360 */
    233
    234        /**
    235         * Pre-Calculates the values (used for the adams? series) which will be used to calculate the authalic latitude.
    236         * Snyder (3-18).
    237         *
    238         * @param eccentricitySquared
    239         *            the squared eccentricity from the ellipsoid to calculate the authalic latitude for.
    240         * @return the precalculated values.
    241         * @deprecated use {@link ProjectionUtils#getAuthalicLatitudeSeriesValues(double)} instead.;
    242         */
    243        @Deprecated
    244        public static double[] authset( double eccentricitySquared ) {
    245            return getAuthalicLatitudeSeriesValues( eccentricitySquared );
    246        }
    247
    248        /**
    249         * Pre-Calculates the values (used for the adams? series) which will be used to calculate the authalic latitude.
    250         * Snyder (3-18).
    251         *
    252         * @param eccentricitySquared
    253         *            the squared eccentricity from the ellipsoid to calculate the authalic latitude for.
    254         * @return the precalculated values [0] = e^2/3 + e^4*(31/180) + e^6*(517/5040), [1]= e^4*(23/360) + e^6*(251/3780)
    255         *         and [2] = e^6*(761/45360).
    256         */
    257        public static double[] getAuthalicLatitudeSeriesValues( double eccentricitySquared ) {
    258            double[] precalculatedSerie = new double[3];
    259            precalculatedSerie[0] = eccentricitySquared * P00;
    260            double t = eccentricitySquared * eccentricitySquared;
    261            precalculatedSerie[0] += t * P01;
    262            precalculatedSerie[1] = t * P10;
    263            t *= eccentricitySquared;
    264            precalculatedSerie[0] += t * P02;
    265            precalculatedSerie[1] += t * P11;
    266            precalculatedSerie[2] = t * P20;
    267            return precalculatedSerie;
    268        }
    269
    270        /**
    271         * Gets phi from the authalic latitude beta and the precalculated values of the adams? serie. From Snyder (3-18).
    272         *
    273         * @param beta
    274         *            authalic latitude.
    275         * @param APA
    276         *            the precalculated values from the series gotten from
    277         *            {@link ProjectionUtils#getAuthalicLatitudeSeriesValues(double)}.
    278         * @return the phi on the ellipsoid.
    279         * @deprecated use {@link ProjectionUtils#calcPhiFromAuthalicLatitude(double, double[])} instead.;
    280         */
    281        @Deprecated
    282        public static double authlat( double beta, double[] APA ) {
    283            return calcPhiFromAuthalicLatitude( beta, APA );
    284        }
    285
    286        /**
    287         * Gets phi from the authalic latitude beta and the precalculated values of the adams? serie. From Snyder (3-18).
    288         *
    289         * @param beta
    290         *            authalic latitude.
    291         * @param APA
    292         *            the precalculated values from the serie gotten from
    293         *            {@link ProjectionUtils#getAuthalicLatitudeSeriesValues(double)}.
    294         * @return the phi on the ellipsoid.
    295         */
    296        public static double calcPhiFromAuthalicLatitude( double beta, double[] APA ) {
    297            double t = beta + beta;
    298            return ( beta + APA[0] * Math.sin( t ) + APA[1] * Math.sin( t + t ) + APA[2] * Math.sin( t + t + t ) );
    299        }
    300
    301        /**
    302         * Calcs the length of a vector given by two points x and y
    303         *
    304         * @param dx
    305         *            of the vector
    306         * @param dy
    307         *            of the vector
    308         * @return the length
    309         */
    310        public static double length( double dx, double dy ) {
    311            return Math.hypot( dx, dy );
    312        }
    313
    314        /**
    315         * This method calculates the innerpart of the conformal latitude's definition (Snyder p.15 3-1). This formula is
    316         * almost equal to the calculation of the half colatitude from the conformal latitude (Snyder p.108 15-9). They only
    317         * differ a sign in the first term.
    318         *
    319         * @param phi
    320         *            to calculate the conformal latitude from
    321         * @param sinphi
    322         *            the sinus of the phi.
    323         * @param eccentricity
    324         *            of the ellipsoid to which the phi should be made conformal to.
    325         * @return the value of the innerpart of the conformal latitude formula. i.e. tan( pi/4 <b>+</b> phi/2)<b>*</b>[(1-e*sin(phi))/1+e*sin(phi))]^e/2.
    326         */
    327        public static double conformalLatitudeInnerPart( double phi, double sinphi, double eccentricity ) {
    328            sinphi *= eccentricity;
    329            return ( Math.tan( .5 * ( HALFPI + phi ) ) ) * Math.pow( ( 1. - sinphi ) / ( 1. + sinphi ), .5 * eccentricity );
    330        }
    331
    332        /**
    333         * This method calculates the innerpart of the conformal latitude's definition (Snyder p.15 3-1). This formula is
    334         * almost equal to the calculation of the half colatitude from the conformal latitude (Snyder p.108 15-9). They only
    335         * differ a sign in the first term.
    336         *
    337         * @param phi
    338         *            to calculate the conformal latitude from
    339         * @param sinphi
    340         *            the sinus of the phi.
    341         * @param eccentricity
    342         *            of the ellipsoid to which the phi should be made conformal to.
    343         * @return the value of the innerpart of the conformal latitude formula. i.e. tan( pi/4 <b>+</b>
    344         *         phi/2)*[(1-e*sin(phi))/1+e*sin(phi))]^e/2.
    345         * @deprecated Use {@link #conformalLatitudeInnerPart(double,double,double)} instead
    346         */
    347        @Deprecated
    348        public static double ssfn( double phi, double sinphi, double eccentricity ) {
    349            return conformalLatitudeInnerPart( phi, sinphi, eccentricity );
    350        }
    351
    352        /**
    353         * This method calculates the tangens of the half colatitude from the conformal latitude (Snyder p.108 15-9).
    354         *
    355         * @param phi
    356         *            to calculate the half of the co latitude of the conformal latitude from
    357         * @param sinphi
    358         *            the sinus of the phi.
    359         * @param eccentricity
    360         *            of the ellipsoid to which the phi should be made conformal to.
    361         * @return the value of the tangens of half of the conformal latitude formula. i.e. tan( pi/4 <b>-</b> phi/2)<b>/</b>[(1-e*sin(phi))/1+e*sin(phi))]^e/2.
    362         */
    363        public static double tanHalfCoLatitude( double phi, double sinphi, double eccentricity ) {
    364            sinphi *= eccentricity;
    365            return ( Math.tan( .5 * ( HALFPI - phi ) ) ) / Math.pow( ( 1. - sinphi ) / ( 1. + sinphi ), .5 * eccentricity );
    366        }
    367
    368        /**
    369         * This method calculates the tangens of the half colatitude from the conformal latitude (Snyder p.108 15-9). This
    370         * formula is almost equal to the calculation of the innerpart of the conformal latitude's definition (Snyder p.15
    371         * 3-1). They only differ a sign in the first term.
    372         *
    373         * @param phi
    374         *            to calculate the half of the co latitude of the conformal latitude from
    375         * @param sinphi
    376         *            the sinus of the phi.
    377         * @param eccentricity
    378         *            of the ellipsoid to which the phi should be made conformal to.
    379         * @return the value of the innerpart of the conformal latitude formula (given sign + or -). i.e. tan( pi/4 (+-)
    380         *         phi/2)*[(1-e*sin(phi))/1+e*sin(phi))]^e/2.
    381         * @deprecated Use {@link #tanHalfCoLatitude(double,double,double)} instead
    382         */
    383        @Deprecated
    384        public static double tsfn( double phi, double sinphi, double eccentricity ) {
    385            return tanHalfCoLatitude( phi, sinphi, eccentricity );
    386        }
    387
    388        /**
    389         * This method can be used to calculate the value of a variable called 'm' by Snyder (Snyder p.101 14-15).
    390         *
    391         * @param sinphi
    392         *            the sinus of the phi
    393         * @param cosphi
    394         *            the cosinus of the phi
    395         * @param eccentricitySquared
    396         *            the value eccentricity * eccentricity.
    397         * @return cos( phi) / Math.sqrt( 1 - eccentricity*eccentricity*sin(phi)*sin(phi) ).
    398         * @deprecated Use {@link #calcMFromSnyder(double,double,double)} instead
    399         */
    400        @Deprecated
    401        public static double msfn( double sinphi, double cosphi, double eccentricitySquared ) {
    402            return calcMFromSnyder( sinphi, cosphi, eccentricitySquared );
    403        }
    404
    405        /**
    406         * This method can be used to calculate the value of a variable called 'm' by Snyder (Snyder p.101 14-15).
    407         *
    408         * @param sinphi
    409         *            the sinus of the phi
    410         * @param cosphi
    411         *            the cosinus of the phi
    412         * @param eccentricitySquared
    413         *            the value eccentricity * eccentricity.
    414         * @return cos( phi) / Math.sqrt( 1 - eccentricity*eccentricity*sin(phi)*sin(phi) ).
    415         */
    416        public static double calcMFromSnyder( double sinphi, double cosphi, double eccentricitySquared ) {
    417            return cosphi / Math.sqrt( 1.0 - eccentricitySquared * sinphi * sinphi );
    418        }
    419
    420        /**
    421         * Copied these value from proj4, I think they are reformed to fit some rule... but I don't know which rule that is
    422         * :-(
    423         */
    424        private final static double C00 = 1; /* 1 :-) */
    425
    426        private final static double C02 = .25; /* 1/4 */
    427
    428        private final static double C04 = .046875;/* 3/64 */
    429
    430        private final static double C06 = .01953125;/* 5/256 */
    431
    432        private final static double C08 = .01068115234375;/* 175 / 16384 */
    433
    434        private final static double C22 = .75;
    435
    436        private final static double C44 = .46875;
    437
    438        private final static double C46 = .01302083333333333333;
    439
    440        private final static double C48 = .00712076822916666666;
    441
    442        private final static double C66 = .36458333333333333333;
    443
    444        private final static double C68 = .00569661458333333333;
    445
    446        private final static double C88 = .3076171875;
    447
    448        /**
    449         * Pre Calculates the values for the series to calculate for a given ellipsoid with given eccentricity the distance
    450         * along the meridian from the equator to a given latitude
    451         * {@link #getDistanceAlongMeridian(double, double, double, double[])}.
    452         *
    453         * @param es
    454         *            the squared eccentricity of the underlying ellipsoid.
    455         * @return the precalculated values for given ellipsoid.
    456         * @deprecated Use {@link #getRectifiyingLatitudeValues(double)} instead
    457         */
    458        @Deprecated
    459        public static double[] enfn( double es ) {
    460            return getRectifiyingLatitudeValues( es );
    461        }
    462
    463        /**
    464         * Pre Calculates the values for the series to calculate for a given ellipsoid with given eccentricity the distance
    465         * along the meridian from the equator to a given latitude
    466         * {@link #getDistanceAlongMeridian(double, double, double, double[])}.
    467         *
    468         * @param es
    469         *            the squared eccentricity of the underlying ellipsoid.
    470         * @return the precalculated values for given ellipsoid.
    471         */
    472        public static double[] getRectifiyingLatitudeValues( double es ) {
    473            double[] en = new double[5];
    474            en[0] = C00 - es * ( C02 + es * ( C04 + es * ( C06 + es * C08 ) ) );
    475            en[1] = es * ( C22 - es * ( C04 + es * ( C06 + es * C08 ) ) );
    476
    477            double t = es * es;
    478            en[2] = t * ( C44 - es * ( C46 + es * C48 ) );
    479
    480            t *= es;
    481            en[3] = t * ( C66 - es * C68 );
    482
    483            en[4] = t * es * C88;
    484
    485            return en;
    486        }
    487
    488        /**
    489         * This method calcs for a a given ellispoid the distance along the meridian from the equator to latitude phi Snyder
    490         * (p.17 3-21). It is used to calculate the rectifying latitude <i>mu</i>.
    491         *
    492         * @param phi
    493         *            the lattitude of the point in radians
    494         * @param sphi
    495         *            the sinus of the latitude
    496         * @param cphi
    497         *            the cosinus of the latitude
    498         * @param en
    499         *            an array (of length 5) containing the precalculate values for this ellipsoid gotten from
    500         *            {@link #getRectifiyingLatitudeValues(double)}.
    501         * @return the distance along the meridian from the equator to latitude phi.
    502         * @deprecated Use {@link #getDistanceAlongMeridian(double,double,double,double[])} instead
    503         */
    504        @Deprecated
    505        public static double mlfn( double phi, double sphi, double cphi, double[] en ) {
    506            return getDistanceAlongMeridian( phi, sphi, cphi, en );
    507        }
    508
    509        /**
    510         * This method calcs the distance along the meridian from the equator to latitude phi for a a given ellispoid Snyder
    511         * (p.17 3-21). It is used to calculate the rectifying latitude <i>mu</i>.
    512         *
    513         * @param phi
    514         *            the lattitude of the point in radians
    515         * @param sphi
    516         *            the sinus of the latitude
    517         * @param cphi
    518         *            the cosinus of the latitude
    519         * @param en
    520         *            an array (of length 5) containing the precalculate values for this ellipsoid gotten from
    521         *            {@link #getRectifiyingLatitudeValues(double)}.
    522         * @return the distance along the meridian from the equator to latitude phi.
    523         */
    524        public static double getDistanceAlongMeridian( double phi, double sphi, double cphi, double[] en ) {
    525            cphi *= sphi;
    526            sphi *= sphi;
    527            return ( en[0] * phi ) - ( cphi * ( en[1] + sphi * ( en[2] + sphi * ( en[3] + sphi * en[4] ) ) ) );
    528        }
    529
    530        /**
    531         * This method calcs lattitude phi from a given distance along the meridian to the equator for a a given ellispoid
    532         * and is therefore the inverse of the {@link #getDistanceAlongMeridian(double, double, double, double[])}. Phi is
    533         * determined to EPS (1e-11) radians, which is about 1e-6 seconds.
    534         *
    535         * @param initialValue
    536         *            to calculate phi from, a good starting value is using the (distance along the meridian / y*scale) e.g.
    537         *            the scaled y value on the meridian.
    538         * @param squaredEccentricity
    539         *            the squared eccentricity of the ellipsoid.
    540         * @param en
    541         *            an array (of length 5) containing the precalculate values for this ellipsoid gotten from
    542         *            {@link #getRectifiyingLatitudeValues(double)}.
    543         * @return the lattitude phi.
    544         * @deprecated Use {@link #calcPhiFromMeridianDistance(double,double,double[])} instead
    545         */
    546        @Deprecated
    547        public static double inv_mlfn( double initialValue, double squaredEccentricity, double[] en ) {
    548            return calcPhiFromMeridianDistance( initialValue, squaredEccentricity, en );
    549        }
    550
    551        /**
    552         * This method calcs lattitude phi from a given distance along the meridian to the equator for a a given ellispoid
    553         * and is therefore the inverse of the {@link #getDistanceAlongMeridian(double, double, double, double[])}. Phi is
    554         * determined to EPS (1e-11) radians, which is about 1e-6 seconds.
    555         *
    556         * @param initialValue
    557         *            to calculate phi from, a good starting value is using the (distance along the meridian / y*scale) e.g.
    558         *            the scaled y value on the meridian.
    559         * @param squaredEccentricity
    560         *            the squared eccentricity of the ellipsoid.
    561         * @param en
    562         *            an array (of length 5) containing the precalculate values for this ellipsoid gotten from
    563         *            {@link #getRectifiyingLatitudeValues(double)}.
    564         * @return the lattitude phi or the best approximated value if no suitable convergence was found.
    565         */
    566        public static double calcPhiFromMeridianDistance( double initialValue, double squaredEccentricity, double[] en ) {
    567            double k = 1. / ( 1. - squaredEccentricity );
    568            double phi = initialValue;
    569            /* (from proj4: ->) rarely goes over 2 iterations */
    570            for ( int i = MAX_ITER; i != 0; i-- ) {
    571                double s = Math.sin( phi );
    572                double t = 1. - squaredEccentricity * s * s;
    573                t = ( getDistanceAlongMeridian( phi, s, Math.cos( phi ), en ) - initialValue ) * ( t * Math.sqrt( t ) ) * k;
    574                phi -= t;
    575                if ( Math.abs( t ) < EPS11 ) {
    576                    return phi;
    577                }
    578            }
    579            return phi;
    580        }
    581
    582        /**
    583         * A helper method, which returns the acos from value or if value < -1 pi or value>1 0.
    584         *
    585         * @param value
    586         *            (in radians) from which the acos must be calculated
    587         * @return the acos from value or if value < -1 pi or if value > 1 0.
    588         */
    589        public static double acosScaled( double value ) {
    590            return ( value < -1 ) ? Math.PI : ( value > 1 ) ? 0 : Math.acos( value );
    591        }
    592
    593        /**
    594         * A helper method, which returns the asin from value or if value < -1 (-pi/2) or value>1 (pi/2).
    595         *
    596         * @param value
    597         *            (in radians) from which the asin must be calculated
    598         * @return the asin from value or if value < -1 (-pi/2) or value>1 (pi/2).
    599         */
    600        public static double asinScaled( double value ) {
    601            return ( value < -1 ) ? -HALFPI : ( value > 1 ) ? HALFPI : Math.asin( value );
    602        }
    603
    604        /**
    605         * A helper method modulos (pi)the given angle (in radians) until the result fits betwee -HALFPI and HALF_PI.
    606         *
    607         * @param angle
    608         *            in radians
    609         * @return the angle adjusted to -pi/2 + pi/2 or 0 if the angle is NaN or Infinite.
    610         */
    611        public static double normalizeLatitude( double angle ) {
    612            if ( Double.isInfinite( angle ) || Double.isNaN( angle ) ) {
    613                return 0;
    614            }
    615            while ( angle > HALFPI ) {
    616                angle -= Math.PI;
    617            }
    618            while ( angle < -HALFPI ) {
    619                angle += Math.PI;
    620            }
    621            return angle;
    622        }
    623
    624        /**
    625         * A helper method modulos (2*pi)the given angle (in radians) until the result fits betwee -PI and PI.
    626         *
    627         * @param angle
    628         *            to be normalized
    629         * @return the angle adjusted to -2*pi + pi*2 or 0 if the angle is NaN or Infinite.
    630         */
    631        public static double normalizeLongitude( double angle ) {
    632            if ( Double.isInfinite( angle ) || Double.isNaN( angle ) ) {
    633                return 0;
    634            }
    635            while ( angle > Math.PI ) {
    636                angle -= TWOPI;
    637            }
    638            while ( angle < -Math.PI ) {
    639                angle += TWOPI;
    640            }
    641            return angle;
    642        }
    643
    644        /**
    645         * Converts a Deegree.MinSec value into it's radian equivalent. <code>
    646         * for example 13.120637 dms -> 13.201769444444446° -> 0.23041434389473822 rd
    647         * </code>
    648         *
    649         * @param inCoord
    650         *            to be converted to radians.
    651         * @return the radian equivalent of the inCoord.
    652         */
    653        public static double DecMinSecToRadians( double inCoord ) {
    654            // get decimal minutes
    655            double remainder = getRemainder( inCoord ) * 100;
    656            // add the decimal minutes to the decimal seconds
    657            remainder = ( Math.floor( remainder ) * 60 ) + ( getRemainder( remainder ) * 100 );
    658            return ( Math.floor( inCoord ) + ( remainder / 3600 ) ) * DTR;
    659        }
    660
    661        /**
    662         * Converts a radian to its Deegree.MinSec equivalent.<code>
    663         * For example 0.23041434389473822 rd -> 13.201769444444446° -> 13.120637 dms
    664         * </code>
    665         *
    666         * @param inCoord
    667         *            to be converted to degrees.minsec
    668         * @return the radian equivalent of the inCoord.
    669         */
    670        public static double RadiansToDecMinSec( double inCoord ) {
    671
    672            double degrees = RTD * inCoord;
    673            double remainder = getRemainder( degrees );
    674            return Math.floor( degrees ) + ( Math.floor( remainder * 60 ) / 100 )
    675                   + ( ( remainder * .36 ) - ( Math.floor( remainder * 60 ) * .006 ) );
    676        }
    677
    678        /**
    679         * Retrieve the remainder of a given double value, e.g. value - Math.floor( value ).
    680         *
    681         * @param value
    682         *            to get the remainder from.
    683         *
    684         * @return the remainder of the given value.
    685         */
    686        private static double getRemainder( double value ) {
    687            return value - Math.floor( value );
    688        }
    689    }