thoerup/datumconversion/DatumConverter.java

package dk.thoerup.datumconversion;

public class DatumConverter {

        static Utm LLtoUTM(int ReferenceEllipsoid, Ll ll)
        {
                //converts lat/long to UTM coords.  Equations from USGS Bulletin 1532 
                //East Longitudes are positive, West longitudes are negative. 
                //North latitudes are positive, South latitudes are negative
                //Lat and Long are in decimal degrees
                //Written by Chuck Gantz- chuck.gantz@globalstar.com


                double UTMNorthing; 
                double UTMEasting; 
                String UTMZone;
                
                double Lat = ll.lattitude;
                double Long = ll.longitude;


                double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius;
                double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared;
                double k0 = 0.9996;

                double LongOrigin;
                double eccPrimeSquared;
                double N, T, C, A, M;

                //Make sure the longitude is between -180.00 .. 179.9
                double LongTemp = (Long+180) - ((int)(Long+180)/360)*360-180; // -180.00 .. 179.9;

                double LatRad = Lat * Constants.deg2rad;
                double LongRad = LongTemp * Constants.deg2rad;
                double LongOriginRad;
                int    ZoneNumber;

                ZoneNumber = ( (int)((LongTemp + 180)/6) ) + 1;

                if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
                        ZoneNumber = 32;

                // Special zones for Svalbard
                if( Lat >= 72.0 && Lat < 84.0 ) 
                {
                        if(      LongTemp >= 0.0  && LongTemp <  9.0 ) ZoneNumber = 31;
                        else if( LongTemp >= 9.0  && LongTemp < 21.0 ) ZoneNumber = 33;
                        else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35;
                        else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37;
                }
                LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
                LongOriginRad = LongOrigin * Constants.deg2rad;

                //compute the UTM Zone from the latitude and longitude
                UTMZone = ""+ ZoneNumber + UTMLetterDesignator(Lat);

                eccPrimeSquared = (eccSquared)/(1-eccSquared);

                N = a/Math.sqrt(1-eccSquared * Math.sin(LatRad) * Math.sin(LatRad));
                T = Math.tan(LatRad) * Math.tan(LatRad);
                C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad);
                A = Math.cos(LatRad)*(LongRad-LongOriginRad);

                M = a*((1       - eccSquared/4          - 3*eccSquared*eccSquared/64    - 5*eccSquared*eccSquared*eccSquared/256)*LatRad 
                                - (3*eccSquared/8       + 3*eccSquared*eccSquared/32    + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatRad)
                                + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatRad) 
                                - (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatRad));

                UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
                                + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
                                + 500000.0);

                UTMNorthing = (double)(k0*(M+N*Math.tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
                                + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
                if(Lat < 0)
                        UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere

                return new Utm(UTMNorthing, UTMEasting,UTMZone);
        }

        static char UTMLetterDesignator(double Lat)
        {
                //This routine determines the correct UTM letter designator for the given latitude
                //returns 'Z' if latitude is outside the UTM limits of 84N to 80S
                //Written by Chuck Gantz- chuck.gantz@globalstar.com
                char LetterDesignator;

                if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X';
                else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W';
                else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V';
                else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U';
                else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T';
                else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S';
                else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R';
                else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q';
                else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P';
                else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N';
                else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M';
                else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L';
                else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K';
                else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J';
                else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H';
                else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G';
                else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F';
                else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E';
                else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D';
                else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C';
                else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits

                return LetterDesignator;
        }


        static Ll UTMtoLL(int ReferenceEllipsoid, Utm utm )
        {
                //converts UTM coords to lat/long.  Equations from USGS Bulletin 1532 
                //East Longitudes are positive, West longitudes are negative. 
                //North latitudes are positive, South latitudes are negative
                //Lat and Long are in decimal degrees. 
                //Written by Chuck Gantz- chuck.gantz@globalstar.com

                double k0 = 0.9996;
                double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius;
                double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared;
                double eccPrimeSquared;
                double e1 = (1-Math.sqrt(1-eccSquared))/(1+Math.sqrt(1-eccSquared));
                double N1, T1, C1, R1, D, M;
                double LongOrigin;
                double mu, phi1, phi1Rad;
                double x, y;
                int ZoneNumber;
                char ZoneLetter;
                int NorthernHemisphere; //1 for northern hemispher, 0 for southern

                double Lat;
                double Long;
                
                double UTMNorthing = utm.northing;
                double UTMEasting = utm.easting;
                String UTMZone = utm.utmZone;
                
                x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
                y = UTMNorthing;

                ZoneNumber = Integer.parseInt(UTMZone.substring(0,UTMZone.length()-1));
                ZoneLetter = UTMZone.charAt(UTMZone.length()-1);
                //ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
                if((ZoneLetter - 'N') >= 0)
                        NorthernHemisphere = 1;//point is in northern hemisphere
                else
                {
                        NorthernHemisphere = 0;//point is in southern hemisphere
                        y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
                }

                LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone

                eccPrimeSquared = (eccSquared)/(1-eccSquared);

                M = y / k0;
                mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));

                phi1Rad = mu    + (3*e1/2-27*e1*e1*e1/32)*Math.sin(2*mu) 
                + (21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(4*mu)
                +(151*e1*e1*e1/96)*Math.sin(6*mu);
                phi1 = phi1Rad*Constants.rad2deg;

                N1 = a/Math.sqrt(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad));
                T1 = Math.tan(phi1Rad)*Math.tan(phi1Rad);
                C1 = eccPrimeSquared*Math.cos(phi1Rad)*Math.cos(phi1Rad);
                R1 = a*(1-eccSquared)/Math.pow(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad), 1.5);
                D = x/(N1*k0);

                Lat = phi1Rad - (N1*Math.tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
                                +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
                Lat = Lat * Constants.rad2deg;

                Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
                                *D*D*D*D*D/120)/Math.cos(phi1Rad);
                Long = LongOrigin + Long * Constants.rad2deg;

                return new Ll(Lat,Long);
        }

        /////////////////////////////////////////////////////////////////////////////////////////////////////// 
        
        static Swiss LLtoSwissGrid(double Lat, double Long)
        {
        //converts lat/long to Swiss Grid coords.  Equations from "Supplementary PROJ.4 Notes-
        //Swiss Oblique Mercator Projection", August 5, 1995, Release 4.3.3, by Gerald I. Evenden
        //Lat and Long are in decimal degrees
        //This transformation is, of course, only valid in Switzerland
                //Written by Chuck Gantz- chuck.gantz@globalstar.com
                double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid
                double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
                double ecc = Math.sqrt(eccSquared);

                double LongOrigin = 7.43958333; //E7d26'22.500"
                double LatOrigin = 46.95240556; //N46d57'8.660"

                double SwissNorthing;
                double SwissEasting;
                
                double LatRad = Lat*Constants.deg2rad;
                double LongRad = Long*Constants.deg2rad;
                double LatOriginRad = LatOrigin*Constants.deg2rad;
                double LongOriginRad = LongOrigin*Constants.deg2rad;

                double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); 
                
                double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);

                //eqn. 1
                double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2))
                                        -c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) 
                                                - ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad))));

                
                double LongRadPrime = c*(LongRad - LongOriginRad); //eqn 2
                double w = c*(Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) 
                                        - ecc/2 * Math.log((1+ecc*Math.sin(LatRad)) / (1-ecc*Math.sin(LatRad)))) + K; //eqn 1
                double LatRadPrime = 2 * (Math.atan(Math.exp(w)) - Constants.FOURTHPI); //eqn 1
                
                //eqn 3
                double sinLatDoublePrime = Math.cos(equivLatOrgRadPrime) * Math.sin(LatRadPrime) 
                                                                        - Math.sin(equivLatOrgRadPrime) * Math.cos(LatRadPrime) * Math.cos(LongRadPrime);
                double LatRadDoublePrime = Math.asin(sinLatDoublePrime);
                
                //eqn 4
                double sinLongDoublePrime = Math.cos(LatRadPrime)*Math.sin(LongRadPrime) / Math.cos(LatRadDoublePrime);
                double LongRadDoublePrime = Math.asin(sinLongDoublePrime);

                double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad));

                SwissNorthing = R*Math.log(Math.tan(Constants.FOURTHPI + LatRadDoublePrime/2)) + 200000.0; //eqn 5
                SwissEasting = R*LongRadDoublePrime + 600000.0; //eqn 6

                return new Swiss(SwissNorthing, SwissEasting);
        }


        static Ll SwissGridtoLL(double SwissNorthing, double SwissEasting)
        {
                double Lat;
                double Long;
                
                double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid
                double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
                double ecc = Math.sqrt(eccSquared);

                double LongOrigin = 7.43958333; //E7d26'22.500"
                double LatOrigin = 46.95240556; //N46d57'8.660"

                double LatOriginRad = LatOrigin*Constants.deg2rad;
                double LongOriginRad = LongOrigin*Constants.deg2rad;

                double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad));

                double LatRadDoublePrime = 2*(Math.atan(Math.exp((SwissNorthing - 200000.0)/R)) - Constants.FOURTHPI); //eqn. 7
                double LongRadDoublePrime = (SwissEasting - 600000.0)/R; //eqn. 8 with equation corrected


                double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); 
                double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);

                double sinLatRadPrime = Math.cos(equivLatOrgRadPrime)*Math.sin(LatRadDoublePrime)
                                                                + Math.sin(equivLatOrgRadPrime)*Math.cos(LatRadDoublePrime)*Math.cos(LongRadDoublePrime);
                double LatRadPrime = Math.asin(sinLatRadPrime);

                double sinLongRadPrime = Math.cos(LatRadDoublePrime)*Math.sin(LongRadDoublePrime)/Math.cos(LatRadPrime);
                double LongRadPrime = Math.asin(sinLongRadPrime);

                Long = (LongRadPrime/c + LongOriginRad) * Constants.rad2deg;

                Lat = NewtonRaphson(LatRadPrime) * Constants.rad2deg;

                return new Ll(Lat, Long);
        }

        static double NewtonRaphson( double initEstimate)
        {
                double Estimate = initEstimate;
                double tol = 0.00001;
                double corr;

                double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
                double ecc = Math.sqrt(eccSquared);

                double LatOrigin = 46.95240556; //N46d57'8.660"
                double LatOriginRad = LatOrigin*Constants.deg2rad;

                double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); 
                
                double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);

                //eqn. 1
                double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2))
                                        -c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) 
                                                - ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad))));
                double C = (K - Math.log(Math.tan(Constants.FOURTHPI + initEstimate/2)))/c;

                do
                {
                        corr = CorrRatio(Estimate, C);
                        Estimate = Estimate - corr;
                }
                while (Math.abs(corr) > tol);

                return Estimate;
        }


        static double CorrRatio(double LatRad, double C)
        {
                double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
                double ecc = Math.sqrt(eccSquared);
                double corr = (C + Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) 
                                        - ecc/2 * Math.log((1+ecc*Math.sin(LatRad)) / (1-ecc*Math.sin(LatRad)))) * (((1-eccSquared*Math.sin(LatRad)*Math.sin(LatRad)) * Math.cos(LatRad)) / (1-eccSquared));

                return corr;
        }
                
        
}
1	torben	262	package dk.thoerup.datumconversion;
2
3			public class DatumConverter {
4
5			static Utm LLtoUTM(int ReferenceEllipsoid, Ll ll)
6			{
7			//converts lat/long to UTM coords. Equations from USGS Bulletin 1532
8			//East Longitudes are positive, West longitudes are negative.
9			//North latitudes are positive, South latitudes are negative
10			//Lat and Long are in decimal degrees
11			//Written by Chuck Gantz- chuck.gantz@globalstar.com
12
13
14			double UTMNorthing;
15			double UTMEasting;
16			String UTMZone;
17
18			double Lat = ll.lattitude;
19			double Long = ll.longitude;
20
21
22
23			double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius;
24			double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared;
25			double k0 = 0.9996;
26
27			double LongOrigin;
28			double eccPrimeSquared;
29			double N, T, C, A, M;
30
31			//Make sure the longitude is between -180.00 .. 179.9
32			double LongTemp = (Long+180) - ((int)(Long+180)/360)*360-180; // -180.00 .. 179.9;
33
34			double LatRad = Lat * Constants.deg2rad;
35			double LongRad = LongTemp * Constants.deg2rad;
36			double LongOriginRad;
37			int ZoneNumber;
38
39			ZoneNumber = ( (int)((LongTemp + 180)/6) ) + 1;
40
41			if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
42			ZoneNumber = 32;
43
44			// Special zones for Svalbard
45			if( Lat >= 72.0 && Lat < 84.0 )
46			{
47			if( LongTemp >= 0.0 && LongTemp < 9.0 ) ZoneNumber = 31;
48			else if( LongTemp >= 9.0 && LongTemp < 21.0 ) ZoneNumber = 33;
49			else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35;
50			else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37;
51			}
52			LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
53			LongOriginRad = LongOrigin * Constants.deg2rad;
54
55			//compute the UTM Zone from the latitude and longitude
56			UTMZone = ""+ ZoneNumber + UTMLetterDesignator(Lat);
57
58			eccPrimeSquared = (eccSquared)/(1-eccSquared);
59
60			N = a/Math.sqrt(1-eccSquared * Math.sin(LatRad) * Math.sin(LatRad));
61			T = Math.tan(LatRad) * Math.tan(LatRad);
62			C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad);
63			A = Math.cos(LatRad)*(LongRad-LongOriginRad);
64
65			M = a((1 - eccSquared/4 - 3eccSquaredeccSquared/64 - 5eccSquaredeccSquaredeccSquared/256)*LatRad
66			- (3eccSquared/8 + 3eccSquaredeccSquared/32 + 45eccSquaredeccSquaredeccSquared/1024)Math.sin(2LatRad)
67			+ (15eccSquaredeccSquared/256 + 45eccSquaredeccSquaredeccSquared/1024)Math.sin(4*LatRad)
68			- (35eccSquaredeccSquaredeccSquared/3072)Math.sin(6*LatRad));
69
70			UTMEasting = (double)(k0N(A+(1-T+C)AA*A/6
71			+ (5-18T+TT+72C-58eccPrimeSquared)AAAA*A/120)
72			+ 500000.0);
73
74			UTMNorthing = (double)(k0(M+NMath.tan(LatRad)(AA/2+(5-T+9C+4CC)AAA*A/24
75			+ (61-58T+TT+600C-330eccPrimeSquared)AAAAAA/720)));
76			if(Lat < 0)
77			UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
78
79			return new Utm(UTMNorthing, UTMEasting,UTMZone);
80			}
81
82			static char UTMLetterDesignator(double Lat)
83			{
84			//This routine determines the correct UTM letter designator for the given latitude
85			//returns 'Z' if latitude is outside the UTM limits of 84N to 80S
86			//Written by Chuck Gantz- chuck.gantz@globalstar.com
87			char LetterDesignator;
88
89			if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X';
90			else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W';
91			else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V';
92			else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U';
93			else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T';
94			else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S';
95			else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R';
96			else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q';
97			else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P';
98			else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N';
99			else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M';
100			else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L';
101			else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K';
102			else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J';
103			else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H';
104			else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G';
105			else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F';
106			else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E';
107			else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D';
108			else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C';
109			else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
110
111			return LetterDesignator;
112			}
113
114
115			static Ll UTMtoLL(int ReferenceEllipsoid, Utm utm )
116			{
117			//converts UTM coords to lat/long. Equations from USGS Bulletin 1532
118			//East Longitudes are positive, West longitudes are negative.
119			//North latitudes are positive, South latitudes are negative
120			//Lat and Long are in decimal degrees.
121			//Written by Chuck Gantz- chuck.gantz@globalstar.com
122
123			double k0 = 0.9996;
124			double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius;
125			double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared;
126			double eccPrimeSquared;
127			double e1 = (1-Math.sqrt(1-eccSquared))/(1+Math.sqrt(1-eccSquared));
128			double N1, T1, C1, R1, D, M;
129			double LongOrigin;
130			double mu, phi1, phi1Rad;
131			double x, y;
132			int ZoneNumber;
133			char ZoneLetter;
134			int NorthernHemisphere; //1 for northern hemispher, 0 for southern
135
136			double Lat;
137			double Long;
138
139			double UTMNorthing = utm.northing;
140			double UTMEasting = utm.easting;
141			String UTMZone = utm.utmZone;
142
143			x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
144			y = UTMNorthing;
145
146			ZoneNumber = Integer.parseInt(UTMZone.substring(0,UTMZone.length()-1));
147			ZoneLetter = UTMZone.charAt(UTMZone.length()-1);
148			//ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
149			if((ZoneLetter - 'N') >= 0)
150			NorthernHemisphere = 1;//point is in northern hemisphere
151			else
152			{
153			NorthernHemisphere = 0;//point is in southern hemisphere
154			y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
155			}
156
157			LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
158
159			eccPrimeSquared = (eccSquared)/(1-eccSquared);
160
161			M = y / k0;
162			mu = M/(a(1-eccSquared/4-3eccSquaredeccSquared/64-5eccSquaredeccSquaredeccSquared/256));
163
164			phi1Rad = mu + (3e1/2-27e1e1e1/32)Math.sin(2mu)
165			+ (21e1e1/16-55e1e1e1e1/32)Math.sin(4mu)
166			+(151e1e1e1/96)Math.sin(6*mu);
167			phi1 = phi1Rad*Constants.rad2deg;
168
169			N1 = a/Math.sqrt(1-eccSquaredMath.sin(phi1Rad)Math.sin(phi1Rad));
170			T1 = Math.tan(phi1Rad)*Math.tan(phi1Rad);
171			C1 = eccPrimeSquaredMath.cos(phi1Rad)Math.cos(phi1Rad);
172			R1 = a(1-eccSquared)/Math.pow(1-eccSquaredMath.sin(phi1Rad)*Math.sin(phi1Rad), 1.5);
173			D = x/(N1*k0);
174
175			Lat = phi1Rad - (N1Math.tan(phi1Rad)/R1)(DD/2-(5+3T1+10C1-4C1C1-9eccPrimeSquared)DDDD/24
176			+(61+90T1+298C1+45T1T1-252eccPrimeSquared-3C1C1)DDDDD*D/720);
177			Lat = Lat * Constants.rad2deg;
178
179			Long = (D-(1+2T1+C1)DDD/6+(5-2C1+28T1-3C1C1+8eccPrimeSquared+24T1*T1)
180			DDDD*D/120)/Math.cos(phi1Rad);
181			Long = LongOrigin + Long * Constants.rad2deg;
182
183			return new Ll(Lat,Long);
184			}
185
186			///////////////////////////////////////////////////////////////////////////////////////////////////////
187
188			static Swiss LLtoSwissGrid(double Lat, double Long)
189			{
190			//converts lat/long to Swiss Grid coords. Equations from "Supplementary PROJ.4 Notes-
191			//Swiss Oblique Mercator Projection", August 5, 1995, Release 4.3.3, by Gerald I. Evenden
192			//Lat and Long are in decimal degrees
193			//This transformation is, of course, only valid in Switzerland
194			//Written by Chuck Gantz- chuck.gantz@globalstar.com
195			double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid
196			double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
197			double ecc = Math.sqrt(eccSquared);
198
199			double LongOrigin = 7.43958333; //E7d26'22.500"
200			double LatOrigin = 46.95240556; //N46d57'8.660"
201
202			double SwissNorthing;
203			double SwissEasting;
204
205			double LatRad = Lat*Constants.deg2rad;
206			double LongRad = Long*Constants.deg2rad;
207			double LatOriginRad = LatOrigin*Constants.deg2rad;
208			double LongOriginRad = LongOrigin*Constants.deg2rad;
209
210			double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared)));
211
212			double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);
213
214			//eqn. 1
215			double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2))
216			-c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2))
217			- ecc/2 * Math.log((1+eccMath.sin(LatOriginRad)) / (1-eccMath.sin(LatOriginRad))));
218
219
220			double LongRadPrime = c*(LongRad - LongOriginRad); //eqn 2
221			double w = c*(Math.log(Math.tan(Constants.FOURTHPI + LatRad/2))
222			- ecc/2 * Math.log((1+eccMath.sin(LatRad)) / (1-eccMath.sin(LatRad)))) + K; //eqn 1
223			double LatRadPrime = 2 * (Math.atan(Math.exp(w)) - Constants.FOURTHPI); //eqn 1
224
225			//eqn 3
226			double sinLatDoublePrime = Math.cos(equivLatOrgRadPrime) * Math.sin(LatRadPrime)
227			- Math.sin(equivLatOrgRadPrime) * Math.cos(LatRadPrime) * Math.cos(LongRadPrime);
228			double LatRadDoublePrime = Math.asin(sinLatDoublePrime);
229
230			//eqn 4
231			double sinLongDoublePrime = Math.cos(LatRadPrime)*Math.sin(LongRadPrime) / Math.cos(LatRadDoublePrime);
232			double LongRadDoublePrime = Math.asin(sinLongDoublePrime);
233
234			double R = aMath.sqrt(1-eccSquared) / (1-eccSquaredMath.sin(LatOriginRad) * Math.sin(LatOriginRad));
235
236			SwissNorthing = R*Math.log(Math.tan(Constants.FOURTHPI + LatRadDoublePrime/2)) + 200000.0; //eqn 5
237			SwissEasting = R*LongRadDoublePrime + 600000.0; //eqn 6
238
239			return new Swiss(SwissNorthing, SwissEasting);
240			}
241
242
243			static Ll SwissGridtoLL(double SwissNorthing, double SwissEasting)
244			{
245			double Lat;
246			double Long;
247
248			double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid
249			double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
250			double ecc = Math.sqrt(eccSquared);
251
252			double LongOrigin = 7.43958333; //E7d26'22.500"
253			double LatOrigin = 46.95240556; //N46d57'8.660"
254
255			double LatOriginRad = LatOrigin*Constants.deg2rad;
256			double LongOriginRad = LongOrigin*Constants.deg2rad;
257
258			double R = aMath.sqrt(1-eccSquared) / (1-eccSquaredMath.sin(LatOriginRad) * Math.sin(LatOriginRad));
259
260			double LatRadDoublePrime = 2*(Math.atan(Math.exp((SwissNorthing - 200000.0)/R)) - Constants.FOURTHPI); //eqn. 7
261			double LongRadDoublePrime = (SwissEasting - 600000.0)/R; //eqn. 8 with equation corrected
262
263
264			double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared)));
265			double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);
266
267			double sinLatRadPrime = Math.cos(equivLatOrgRadPrime)*Math.sin(LatRadDoublePrime)
268			+ Math.sin(equivLatOrgRadPrime)Math.cos(LatRadDoublePrime)Math.cos(LongRadDoublePrime);
269			double LatRadPrime = Math.asin(sinLatRadPrime);
270
271			double sinLongRadPrime = Math.cos(LatRadDoublePrime)*Math.sin(LongRadDoublePrime)/Math.cos(LatRadPrime);
272			double LongRadPrime = Math.asin(sinLongRadPrime);
273
274			Long = (LongRadPrime/c + LongOriginRad) * Constants.rad2deg;
275
276			Lat = NewtonRaphson(LatRadPrime) * Constants.rad2deg;
277
278			return new Ll(Lat, Long);
279			}
280
281			static double NewtonRaphson( double initEstimate)
282			{
283			double Estimate = initEstimate;
284			double tol = 0.00001;
285			double corr;
286
287			double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
288			double ecc = Math.sqrt(eccSquared);
289
290			double LatOrigin = 46.95240556; //N46d57'8.660"
291			double LatOriginRad = LatOrigin*Constants.deg2rad;
292
293			double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared)));
294
295			double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c);
296
297			//eqn. 1
298			double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2))
299			-c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2))
300			- ecc/2 * Math.log((1+eccMath.sin(LatOriginRad)) / (1-eccMath.sin(LatOriginRad))));
301			double C = (K - Math.log(Math.tan(Constants.FOURTHPI + initEstimate/2)))/c;
302
303			do
304			{
305			corr = CorrRatio(Estimate, C);
306			Estimate = Estimate - corr;
307			}
308			while (Math.abs(corr) > tol);
309
310			return Estimate;
311			}
312
313
314
315			static double CorrRatio(double LatRad, double C)
316			{
317			double eccSquared = Constants.ellipsoid[3].eccentricitySquared;
318			double ecc = Math.sqrt(eccSquared);
319			double corr = (C + Math.log(Math.tan(Constants.FOURTHPI + LatRad/2))
320			- ecc/2 * Math.log((1+eccMath.sin(LatRad)) / (1-eccMath.sin(LatRad)))) * (((1-eccSquaredMath.sin(LatRad)Math.sin(LatRad)) * Math.cos(LatRad)) / (1-eccSquared));
321
322			return corr;
323			}
324
325
326			}