001: /*
002: * The Unified Mapping Platform (JUMP) is an extensible, interactive GUI
003: * for visualizing and manipulating spatial features with geometry and attributes.
004: *
005: * Copyright (C) 2003 Vivid Solutions
006: *
007: * This program is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU General Public License
009: * as published by the Free Software Foundation; either version 2
010: * of the License, or (at your option) any later version.
011: *
012: * This program is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
015: * GNU General Public License for more details.
016: *
017: * You should have received a copy of the GNU General Public License
018: * along with this program; if not, write to the Free Software
019: * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
020: *
021: * For more information, contact:
022: *
023: * Vivid Solutions
024: * Suite #1A
025: * 2328 Government Street
026: * Victoria BC V8T 5G5
027: * Canada
028: *
029: * (250)385-6040
030: * www.vividsolutions.com
031: */
032:
033: package com.vividsolutions.jump.coordsys;
034:
035: /**
036:
037: * @author $Author: javamap $
038: * @version $Revision: 4 $
039:
040: * <pre>
041: * $Id: Spheroid.java 4 2005-06-16 15:27:48Z javamap $
042: * $Date: 2005-06-16 08:27:48 -0700 (Thu, 16 Jun 2005) $
043: * $Log$
044: * Revision 1.1 2005/06/16 15:24:52 javamap
045: * *** empty log message ***
046: *
047: * Revision 1.2 2005/05/03 15:13:29 javamap
048: * *** empty log message ***
049: *
050: * Revision 1.2 2003/11/05 05:26:52 dkim
051: * Added global header; cleaned up Javadoc.
052: *
053: * Revision 1.1 2003/09/15 20:26:11 jaquino
054: * Reprojection
055: *
056: * Revision 1.2 2003/07/25 17:01:03 gkostadinov
057: * Moved classses reponsible for performing the basic projection to a new
058: * package -- base.
059: *
060: * Revision 1.1 2003/07/24 23:14:43 gkostadinov
061: * adding base projection classes
062: *
063: * Revision 1.1 2003/06/20 18:34:30 gkostadinov
064: * Entering the source code into the CVS.
065: * </pre>
066: *
067: */
068:
069: /**
070: * GRS80 spheroid.
071: */
072: public class Spheroid {
073:
074: public double a;// semimajor axis
075: public double b;// semiminor axis
076: public double f;// flattening
077: public double e;// eccentricity (first)
078:
079: double es;// first eccentricity squared
080: double t1, t2, t3;// t1,...,t6 used for
081: double t4, t5, t6;// area computation (local Math.sinusoidal)
082:
083: public Spheroid(Radius rad) {
084: // constructor base on a and b
085: a = rad.a;
086: if (rad.b > 1.0) {
087: b = rad.b;
088: f = 1.0 - b / a;
089: } else {
090: f = 1.0 / rad.rf;
091: b = a - a * f;
092: }
093: es = f + f - f * f;
094: e = Math.sqrt(es);
095: // constaints for local Math.sinusoidal equal-area projection
096: // used for area computations.
097: double e4;
098: // constaints for local Math.sinusoidal equal-area projection
099: // used for area computations.
100: double e6;
101: // constaints for local Math.sinusoidal equal-area projection
102: // used for area computations.
103: double e8;
104: // constaints for local Math.sinusoidal equal-area projection
105: // used for area computations.
106: double e10;
107: double t0;
108: t0 = a * (1.0 - es);
109: e4 = es * es;
110: e6 = e4 * es;
111: e8 = e6 * es;
112: e10 = e8 * es;
113: t1 = t0
114: * (1.0 + 3.0 * es / 4.0 + 45.0 * e4 / 64.0 + 175.0 * e6
115: / 256.0 + 11025.0 * e8 / 16384.0 + 43659.0 * e10 / 65536.0);
116: t2 = t0
117: * (3.0 * es / 4.0 + 15.0 * e4 / 16.0 + 525.0 * e6
118: / 512.0 + 2205.0 * e8 / 2048.0 + 72765.0 * e10 / 65536.0)
119: / 2.0;
120: t3 = t0
121: * (15.0 * e4 / 64.0 + 105.0 * e6 / 256.0 + 2205.0 * e8
122: / 4096.0 + 10395.0 * e10 / 16384.0) / 4.0;
123: t4 = t0
124: * (35.0 * e6 / 512.0 + 315.0 * e8 / 2048.0 + 31185.0 * e10 / 131072.0)
125: / 6.0;
126: t5 = t0 * (315.0 * e8 / 16384.0 + 3465.0 * e10 / 65536.0) / 8.0;
127: t6 = t0 * (693.0 * e10 / 131072.0) / 10.0;
128: }
129:
130: public double getA() {
131: return a;
132: }
133:
134: public double getB() {
135: return b;
136: }
137:
138: public double getF() {
139: return f;
140: }
141:
142: public double getE() {
143: return e;
144: }
145:
146: public double distance(Geographic r, Geographic s) {
147: // compute Math.sin and cos of latitudes and reduced latitudes
148: double L1;
149: // compute Math.sin and cos of latitudes and reduced latitudes
150: double L2;
151: // compute Math.sin and cos of latitudes and reduced latitudes
152: double sinU1;
153: // compute Math.sin and cos of latitudes and reduced latitudes
154: double sinU2;
155: // compute Math.sin and cos of latitudes and reduced latitudes
156: double cosU1;
157: // compute Math.sin and cos of latitudes and reduced latitudes
158: double cosU2;
159: L1 = Math.atan((1.0 - f) * Math.tan(r.lat));
160: L2 = Math.atan((1.0 - f) * Math.tan(s.lat));
161: sinU1 = Math.sin(L1);
162: sinU2 = Math.sin(L2);
163: cosU1 = Math.cos(L1);
164: cosU2 = Math.cos(L2);
165:
166: // compute delta longitude on the sphere
167: double dl;
168:
169: // compute delta longitude on the sphere
170: double dl1;
171:
172: // compute delta longitude on the sphere
173: double dl2;
174:
175: // compute delta longitude on the sphere
176: double dl3;
177:
178: // compute delta longitude on the sphere
179: double cosdl1;
180:
181: // compute delta longitude on the sphere
182: double sindl1;
183: double cosSigma;
184: double sigma;
185: double azimuthEQ;
186: double tsm;
187: dl = s.lon - r.lon;
188: dl1 = dl;
189: cosdl1 = Math.cos(dl);
190: sindl1 = Math.sin(dl);
191: do {
192: cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
193: sigma = Math.acos(cosSigma);
194: azimuthEQ = Math.asin((cosU1 * cosU2 * sindl1)
195: / Math.sin(sigma));
196: tsm = Math.acos(cosSigma - (2.0 * sinU1 * sinU2)
197: / (Math.cos(azimuthEQ) * Math.cos(azimuthEQ)));
198: dl2 = deltaLongitude(azimuthEQ, sigma, tsm);
199: dl3 = dl1 - (dl + dl2);
200: dl1 = dl + dl2;
201: cosdl1 = Math.cos(dl1);
202: sindl1 = Math.sin(dl1);
203: } while (Math.abs(dl3) > 1.0e-032);
204:
205: // compute expansions A and B
206: double u2;
207:
208: // compute expansions A and B
209: double A;
210:
211: // compute expansions A and B
212: double B;
213: u2 = mu2(azimuthEQ);
214: A = bigA(u2);
215: B = bigB(u2);
216:
217: // compute length of geodesic
218: double dsigma;
219: dsigma = B
220: * Math.sin(sigma)
221: * (Math.cos(tsm) + (B * cosSigma * (-1.0 + 2.0 * (Math
222: .cos(tsm) * Math.cos(tsm)))) / 4.0);
223: return b * (A * (sigma - dsigma));
224: }// END - spheroid::distance
225:
226: public double direction(Geographic r, Geographic s) {
227: // compute Math.sin and cos of latitudes and reduced latitudes
228: double L1;
229: // compute Math.sin and cos of latitudes and reduced latitudes
230: double L2;
231: // compute Math.sin and cos of latitudes and reduced latitudes
232: double sinU1;
233: // compute Math.sin and cos of latitudes and reduced latitudes
234: double sinU2;
235: // compute Math.sin and cos of latitudes and reduced latitudes
236: double cosU1;
237: // compute Math.sin and cos of latitudes and reduced latitudes
238: double cosU2;
239: L1 = Math.atan((1.0 - f) * Math.tan(r.lat));
240: L2 = Math.atan((1.0 - f) * Math.tan(s.lat));
241: sinU1 = Math.sin(L1);
242: sinU2 = Math.sin(L2);
243: cosU1 = Math.cos(L1);
244: cosU2 = Math.cos(L2);
245:
246: // compute delta longitude on the sphere
247: double dl;
248:
249: // compute delta longitude on the sphere
250: double dl1;
251:
252: // compute delta longitude on the sphere
253: double dl2;
254:
255: // compute delta longitude on the sphere
256: double dl3;
257:
258: // compute delta longitude on the sphere
259: double cosdl1;
260:
261: // compute delta longitude on the sphere
262: double sindl1;
263: double cosSigma;
264: double sigma;
265: double azimuthEQ;
266: double tsm;
267: dl = s.lon - r.lon;
268: dl1 = dl;
269: cosdl1 = Math.cos(dl);
270: sindl1 = Math.sin(dl);
271: do {
272: cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
273: sigma = Math.acos(cosSigma);
274: azimuthEQ = Math.asin((cosU1 * cosU2 * sindl1)
275: / Math.sin(sigma));
276: tsm = Math.acos(cosSigma - (2.0 * sinU1 * sinU2)
277: / (Math.cos(azimuthEQ) * Math.cos(azimuthEQ)));
278: dl2 = deltaLongitude(azimuthEQ, sigma, tsm);
279: dl3 = dl1 - (dl + dl2);
280: dl1 = dl + dl2;
281: cosdl1 = Math.cos(dl1);
282: sindl1 = Math.sin(dl1);
283: } while (Math.abs(dl3) > 1.0e-032);
284:
285: // compute expansions A and B
286: double u2;
287:
288: // compute expansions A and B
289: double A;
290:
291: // compute expansions A and B
292: double B;
293: u2 = mu2(azimuthEQ);
294: A = bigA(u2);
295: B = bigB(u2);
296:
297: // compute length of geodesic
298: double dsigma;
299:
300: // compute length of geodesic
301: double d_tmp;
302: dsigma = B
303: * Math.sin(sigma)
304: * (Math.cos(tsm) + (B * cosSigma * (-1.0 + 2.0 * (Math
305: .cos(tsm) * Math.cos(tsm)))) / 4.0);
306: d_tmp = b * (A * (sigma - dsigma));
307:
308: // compute forward azimuth
309: double azimuthFD;
310: azimuthFD = Math.atan2((cosU2 * sindl1), (cosU1 * sinU2 - sinU1
311: * cosU2 * cosdl1));
312: if (azimuthFD < 0.0) {
313: azimuthFD = azimuthFD + 2.0 * Math.PI;
314: }
315: return azimuthFD;
316: }// END - spheroid::direction
317:
318: public Geographic project(Geographic r, double length, double angle) {
319: double e2;
320: double e2s;
321: double L1;
322: double cosU1;
323: double sinU1;
324: double cosa1;
325: double sina1;
326: double sig1;
327: double sinae;
328: double azimuthEQ;
329: double u2;
330: double A;
331: double B;
332: double s1;
333: double sigma;
334: double tsm;
335: double del;
336: double cis;
337: e2 = Math.sqrt(a * a - b * b) / b;// second excentricity
338: e2s = e2 * e2;
339: L1 = Math.atan((1.0 - f) * Math.tan(r.lat));
340: cosU1 = Math.cos(L1);
341: sinU1 = Math.sin(L1);
342: cosa1 = Math.cos(angle);
343: sina1 = Math.sin(angle);
344: sig1 = Math.atan(Math.tan(L1) / cosa1);
345: sinae = cosU1 * sina1;
346: azimuthEQ = Math.asin(sinae);
347: u2 = mu2(azimuthEQ);
348: A = bigA(u2);
349: B = bigB(u2);
350: s1 = length / (b * A);
351: sigma = s1;
352: do {
353: tsm = 2.0 * sig1 + sigma;
354: del = B
355: * Math.sin(sigma)
356: * (Math.cos(tsm) + 0.25
357: * B
358: * Math.cos(sigma)
359: * (-1.0 + 2.0 * (Math.cos(tsm) * Math
360: .cos(tsm))));
361: cis = sigma - (s1 + del);
362: sigma = s1 + del;
363: } while (Math.abs(cis) > 1.0e-032);
364: double cossigma;
365: double sinsigma;
366: cossigma = Math.cos(sigma);
367: sinsigma = Math.sin(sigma);
368:
369: Geographic s = new Geographic();
370: double dm;
371: double dl1;
372: s.lat = sinU1 * cossigma + cosU1 * sinsigma * cosa1;
373: dm = Math.sqrt(sinae * sinae
374: + (sinU1 * sinsigma - cosU1 * cossigma * cosa1)
375: * (sinU1 * sinsigma - cosU1 * cossigma * cosa1));
376: s.lat = Math.atan2(s.lat, ((1.0 - f) * dm));
377: dl1 = Math.atan2((sinsigma * sina1), (cosU1 * cossigma - sinU1
378: * sinsigma * cosa1));
379: s.lon = r.lon + dl1 - deltaLongitude(azimuthEQ, sigma, tsm);
380: return s;
381: }// END - spheroid::project
382:
383: public double meridianRadiusOfCurvature(double latitude) {
384: double er;
385: double el;
386: double M0;
387: er = 1.0 - es * Math.sin(latitude) * Math.sin(latitude);
388: el = Math.pow(er, 1.5);
389: M0 = (a * (1.0 - es)) / el;
390: return M0;
391: }// END - spheroid::meridianRadiusOfCurvature
392:
393: public double primeVerticalRadiusOfCurvature(double latitude) {
394: double T1;
395: double T2;
396: double T3;
397: double N0;
398: T1 = a * a;
399: T2 = T1 * Math.cos(latitude) * Math.cos(latitude);
400: T3 = b * b * Math.sin(latitude) * Math.sin(latitude);
401: N0 = T1 / Math.sqrt(T2 + T3);
402: return N0;
403: }// END - spheroid::primeVerticalRadiusOfCurvature
404:
405: public double deltaLongitude(double azimuth, double sigma,
406: double tsm) {
407: // compute the expansion C
408: double das;
409: // compute the expansion C
410: double C;
411: das = Math.cos(azimuth) * Math.cos(azimuth);
412: C = f / 16.0 * das * (4.0 + f * (4.0 - 3.0 * das));
413: // compute the difference in longitude
414: double ctsm;
415: // compute the difference in longitude
416: double DL;
417: ctsm = Math.cos(tsm);
418: DL = ctsm + C * Math.cos(sigma) * (-1.0 + 2.0 * ctsm * ctsm);
419: DL = sigma + C * Math.sin(sigma) * DL;
420: return (1.0 - C) * f * Math.sin(azimuth) * DL;
421: }// END - spheroid::deltaLongitude
422:
423: public double mu2(double azimuth) {
424: double e2;
425:
426: e2 = Math.sqrt(a * a - b * b) / b;
427: return Math.cos(azimuth) * Math.cos(azimuth) * e2 * e2;
428: }// END - spheroid::mu2
429:
430: public double bigA(double u2) {
431: return 1.0 + u2 / 256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2));
432: }// END - spheroid::bigA
433:
434: public double bigB(double u2) {
435: return u2 / 512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2));
436: }// END - spheroid::bigB
437:
438: public double M(double latitude) {
439: // returns the length of a meridian arc from the equator
440: // to the given latitude (from Richard Rapp).
441: return t1 * latitude - t2 * Math.sin(2.0 * latitude) + t3
442: * Math.sin(4.0 * latitude) - t4
443: * Math.sin(6.0 * latitude) + t5
444: * Math.sin(8.0 * latitude) - t5
445: * Math.sin(10.0 * latitude);
446: }// END - spheroid::M
447:
448: }
|