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