001: /*
002: * <copyright>
003: *
004: * Copyright 2003-2004 BBNT Solutions, LLC
005: * under sponsorship of the Defense Advanced Research Projects
006: * Agency (DARPA).
007: *
008: * You can redistribute this software and/or modify it under the
009: * terms of the Cougaar Open Source License as published on the
010: * Cougaar Open Source Website (www.cougaar.org).
011: *
012: * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
013: * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
014: * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
015: * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
016: * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
017: * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
018: * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
019: * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
020: * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
021: * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
022: * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
023: *
024: * </copyright>
025: */
026: package org.cougaar.logistics.plugin.seanet;
027:
028: /**
029: * A stateful great circle calculator.
030: <p>The theory is presented in
031: <a href="http://stout.bbn.com/~kanderso/java/course/doc/episode-10/index.html">Partial Evaluation - Freeing the essence of a computation</a>.
032: <p>Two lat,lon points are kept in a 3D unit vector representation.
033: <p> Once the two points are set, various distance*() and azimuth() function
034: may be called.
035: <p>The method position() computes the position of a fractional distance along the great circle between the points.
036: */
037: public class GreatCircle {
038:
039: private double a, b, c, d, e, f;
040: private double posLatitude, posLongitude;
041:
042: public GreatCircle() {
043: }
044:
045: public GreatCircle(double lat, double lon) {
046: this .setPoint1(lat, lon);
047: }
048:
049: /** Set the first point.**/
050: public void setPoint1(double lat, double lon) {
051: // Planet2.earth.toV(u("lat"), u("lon"))
052: double ph0 = (0.017453292519943295 * lon);
053: double ph1 = (0.017453292519943295 * lat);
054: double ph2 = Math.tan(ph1);
055: double ph3 = (0.9933056200098587 * ph2);
056: double ph4 = Math.atan(ph3);
057: double ph5 = Math.cos(ph4);
058: double ph6 = Math.cos(ph0);
059: double ph7 = (ph5 * ph6);
060: double ph8 = Math.sin(ph0);
061: double ph9 = (ph5 * ph8);
062: double ph10 = Math.sin(ph4);
063: // return new V(ph7, ph9, ph10);
064: a = ph7;
065: b = ph9;
066: c = ph10;
067: }
068:
069: /** Set the second point.**/
070: public void setPoint2(double lat, double lon) {
071: // Planet2.earth.toV(u("lat"), u("lon"))
072: double ph0 = (0.017453292519943295 * lon);
073: double ph1 = (0.017453292519943295 * lat);
074: double ph2 = Math.tan(ph1);
075: double ph3 = (0.9933056200098587 * ph2);
076: double ph4 = Math.atan(ph3);
077: double ph5 = Math.cos(ph4);
078: double ph6 = Math.cos(ph0);
079: double ph7 = (ph5 * ph6);
080: double ph8 = Math.sin(ph0);
081: double ph9 = (ph5 * ph8);
082: double ph10 = Math.sin(ph4);
083: // return new V(ph7, ph9, ph10);
084: d = ph7;
085: e = ph9;
086: f = ph10;
087: }
088:
089: /** Return the latitude of the current position. **/
090: public double getPositionLatitude() {
091: return this .posLatitude;
092: }
093:
094: /** Return the longiute of the current position. **/
095: public double getPositionLongitude() {
096: return this .posLongitude;
097: }
098:
099: /** Compute distance in nautical miles. **/
100: public double distanceNM(double lat1, double lon1, double lat2,
101: double lon2) {
102: this .setPoint1(lat1, lon1);
103: return distanceNM(lat2, lon2);
104: }
105:
106: /** Distance in nautical miles from the first point and this location. **/
107: public double distanceNM(double lat, double lon) {
108: this .setPoint2(lat, lon);
109: return this .distanceNM();
110: }
111:
112: /** Distance in nautical miles between the two points. **/
113: public double distanceNM() {
114: // http://forum.swarthmore.edu/pow/solutio69.html
115: return distanceKM() / 1.852;
116: }
117:
118: /** Compute distance in kilometers. **/
119: public double distanceKM(double lat, double lon) {
120: this .setPoint2(lat, lon);
121: return this .distanceKM();
122: }
123:
124: /** Distance in kilometers between the two points. **/
125: public double distanceKM() {
126: // times(Planet2.earth.radius(), v1.distance(v2))
127: double ph0 = (e * c);
128: double ph1 = (f * b);
129: double ph2 = (ph0 - ph1);
130: double ph3 = (f * a);
131: double ph4 = (d * c);
132: double ph5 = (ph3 - ph4);
133: double ph6 = (d * b);
134: double ph7 = (e * a);
135: double ph8 = (ph6 - ph7);
136: double ph9 = (ph2 * ph2);
137: double ph10 = (ph5 * ph5);
138: double ph11 = (ph8 * ph8);
139: double ph12 = (ph10 + ph11);
140: double ph13 = (ph9 + ph12);
141: double ph14 = Math.sqrt(ph13);
142: double ph15 = (d * a);
143: double ph16 = (e * b);
144: double ph17 = (f * c);
145: double ph18 = (ph16 + ph17);
146: double ph19 = (ph15 + ph18);
147: double ph20 = Math.atan2(ph14, ph19);
148: double ph21 = (6378.137 * ph20);
149: return ph21;
150: }
151:
152: /** Returns the azimuth from point1 to point2 in degees. **/
153: public double azimuth() {
154: //Planet2.degrees(v1.azimuth(v2))
155:
156: double ph0 = (-b);
157: double ph1 = (ph0 * ph0);
158: double ph2 = (a * a);
159: double ph3 = (ph2 + ph1);
160: double ph4 = Math.sqrt(ph3);
161: boolean ph5 = ph4 < 1.1920929E-7;
162: double ph6;
163: if (ph5) {
164: ph6 = 3.141592653589793;
165: } else {
166: double ph7 = (e * c);
167: double ph8 = (f * b);
168: double ph9 = (ph7 - ph8);
169: double ph10 = (f * a);
170: double ph11 = (d * c);
171: double ph12 = (ph10 - ph11);
172: double ph13 = (d * b);
173: double ph14 = (e * a);
174: double ph15 = (ph13 - ph14);
175: //double ph16 = (- ph15);
176: //double ph17 = (ph0 * ph9);
177: //double ph18 = (a * ph12);
178: //double ph19 = (ph18 + ph17);
179: double ph20 = (-ph15);
180: double ph21 = (ph0 * ph9);
181: double ph22 = (a * ph12);
182: double ph23 = (ph22 + ph21);
183: double ph24 = Math.atan2(ph20, ph23);
184: boolean ph25 = ph24 > 0.0;
185: double ph26;
186: if (ph25) {
187: ph26 = ph24;
188: } else {
189: double ph27 = (6.283185307179586 + ph24);
190: ph26 = ph27;
191: }
192: ph6 = ph26;
193: }
194: double ph28 = (57.29577951308232 * ph6);
195: return ph28;
196: }
197: }
|