gov/nasa/worldwind/globes/EllipsoidalGlobe.java

   1 /*
   2 Copyright (C) 2001, 2006 United States Government
   3 as represented by the Administrator of the
   4 National Aeronautics and Space Administration.
   5 All Rights Reserved.
   6 */
   7 package gov.nasa.worldwind.globes;
   8
   9 import gov.nasa.worldwind.WWObjectImpl;
  10 import gov.nasa.worldwind.geom.*;
  11 import gov.nasa.worldwind.util.Logging;
  12
  13 /**
  14  * @author Tom Gaskins
  15  * @version $Id: EllipsoidalGlobe.java 2570 2007-08-16 22:31:33Z tgaskins $
  16  */
  17 public class EllipsoidalGlobe extends WWObjectImpl implements Globe
  18 {
  19     private final double equatorialRadius;
  20     private final double polarRadius;
  21     private final double es;
  22     private final Vec4 center;
  23
  24     private final ElevationModel elevationModel;
  25
  26     public EllipsoidalGlobe(double equatorialRadius, double polarRadius, double es, ElevationModel em)
  27     {
  28         if (em == null)
  29         {
  30             String msg = Logging.getMessage("nullValue.ElevationModelIsNull");
  31             Logging.logger().severe(msg);
  32             throw new IllegalArgumentException(msg);
  33         }
  34
  35         this.equatorialRadius = equatorialRadius;
  36         this.polarRadius = polarRadius;
  37         this.es = es; // assume it's consistent with the two radii
  38         this.center = Vec4.ZERO;
  39         this.elevationModel = em;
  40     }
  41
  42     public final double getRadius()
  43     {
  44         return this.equatorialRadius;
  45     }
  46
  47     public final double getEquatorialRadius()
  48     {
  49         return this.equatorialRadius;
  50     }
  51
  52     public final double getPolarRadius()
  53     {
  54         return this.polarRadius;
  55     }
  56
  57     public double getMaximumRadius()
  58     {
  59         return this.equatorialRadius;
  60     }
  61
  62     public double getRadiusAt(Angle latitude, Angle longitude)
  63     {
  64         if (latitude == null || longitude == null)
  65         {
  66             String msg = Logging.getMessage("nullValue.AngleIsNull");
  67             Logging.logger().severe(msg);
  68             throw new IllegalArgumentException(msg);
  69         }
  70
  71         return this.computePointFromPosition(latitude, longitude, 0d).getLength3();
  72     }
  73
  74     public double getRadiusAt(LatLon latLon)
  75     {
  76         if (latLon == null)
  77         {
  78             String msg = Logging.getMessage("nullValue.LatLonIsNull");
  79             Logging.logger().severe(msg);
  80             throw new IllegalArgumentException(msg);
  81         }
  82
  83         return this.computePointFromPosition(latLon.getLatitude(), latLon.getLongitude(), 0d).getLength3();
  84     }
  85
  86     public double getEccentricitySquared()
  87     {
  88         return this.es;
  89     }
  90
  91     public final double getDiameter()
  92     {
  93         return this.equatorialRadius * 2;
  94     }
  95
  96     public final Vec4 getCenter()
  97     {
  98         return this.center;
  99     }
 100
 101     public double getMaxElevation()
 102     {
 103         return this.elevationModel.getMaximumElevation();
 104     }
 105
 106     public double getMinElevation()
 107     {
 108         return this.elevationModel.getMinimumElevation();
 109     }
 110
 111     public final Extent getExtent()
 112     {
 113         return this;
 114     }
 115
 116     public boolean intersects(Frustum frustum)
 117     {
 118         if (frustum == null)
 119         {
 120             String message = Logging.getMessage("nullValue.FrustumIsNull");
 121             Logging.logger().severe(message);
 122             throw new IllegalArgumentException(message);
 123         }
 124
 125         return frustum.intersects(this);
 126     }
 127
 128     public Intersection[] intersect(Line line)
 129     {
 130         return this.intersect(line, this.equatorialRadius, this.polarRadius);
 131     }
 132
 133     public Intersection[] intersect(Line line, double altitude)
 134     {
 135         return this.intersect(line, this.equatorialRadius + altitude, this.polarRadius + altitude);
 136     }
 137
 138     private Intersection[] intersect(Line line, double equRadius, double polRadius)
 139     {
 140         if (line == null)
 141         {
 142             String message = Logging.getMessage("nullValue.LineIsNull");
 143             Logging.logger().severe(message);
 144             throw new IllegalArgumentException(message);
 145         }
 146
 147         // Taken from Lengyel, 2Ed., Section 5.2.3, page 148.
 148         double m = equRadius / polRadius;
 149         double n = 1d; //this.equatorialRadius / this.equatorialRadius;
 150         double m2 = m * m;
 151         double n2 = n * n;
 152
 153         double vx = line.getDirection().x;
 154         double vy = line.getDirection().y;
 155         double vz = line.getDirection().z;
 156         double sx = line.getOrigin().x;
 157         double sy = line.getOrigin().y;
 158         double sz = line.getOrigin().z;
 159
 160         double a = vx * vx + m2 * vy * vy + n2 * vz * vz;
 161         double b = 2d * (sx * vx + m2 * sy * vy + n2 * sz * vz);
 162         double c = sx * sx + m2 * sy * sy + n2 * sz * sz - equRadius * polRadius;
 163
 164         double discriminant = discriminant(a, b, c);
 165         if (discriminant < 0)
 166             return null;
 167
 168         double discriminantRoot = Math.sqrt(discriminant);
 169         if (discriminant == 0)
 170         {
 171             Vec4 p = line.getPointAt((-b - discriminantRoot) / (2 * a));
 172             return new Intersection[] {new Intersection(p, true)};
 173         }
 174         else // (discriminant > 0)
 175         {
 176             Vec4 near = line.getPointAt((-b - discriminantRoot) / (2 * a));
 177             Vec4 far = line.getPointAt((-b + discriminantRoot) / (2 * a));
 178             return new Intersection[] {new Intersection(near, false), new Intersection(far, false)};
 179         }
 180     }
 181
 182     static private double discriminant(double a, double b, double c)
 183     {
 184         return b * b - 4 * a * c;
 185     }
 186
 187     public boolean intersects(Line line)
 188     {
 189         if (line == null)
 190         {
 191             String msg = Logging.getMessage("nullValue.LineIsNull");
 192             Logging.logger().severe(msg);
 193             throw new IllegalArgumentException(msg);
 194         }
 195
 196         return line.distanceTo(this.center) <= this.equatorialRadius;
 197     }
 198
 199     public boolean intersects(Plane plane)
 200     {
 201         if (plane == null)
 202         {
 203             String msg = Logging.getMessage("nullValue.PlaneIsNull");
 204             Logging.logger().severe(msg);
 205             throw new IllegalArgumentException(msg);
 206         }
 207
 208         double dq1 = plane.dot(this.center);
 209         return dq1 <= this.equatorialRadius;
 210     }
 211
 212     public Vec4 computeSurfaceNormalAtPoint(Vec4 p)
 213     {
 214         if (p == null)
 215         {
 216             String msg = Logging.getMessage("nullValue.PointIsNull");
 217             Logging.logger().severe(msg);
 218             throw new IllegalArgumentException(msg);
 219         }
 220
 221         p = p.subtract3(this.center);
 222
 223         return new Vec4(
 224             p.x / (this.equatorialRadius * this.equatorialRadius),
 225             p.y / (this.polarRadius * this.polarRadius),
 226             p.z / (this.equatorialRadius * this.equatorialRadius)).normalize3();
 227     }
 228
 229     public final ElevationModel getElevationModel()
 230     {
 231         return this.elevationModel;
 232     }
 233
 234     public final double getElevation(Angle latitude, Angle longitude)
 235     {
 236         if (latitude == null || longitude == null)
 237         {
 238             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 239             Logging.logger().severe(message);
 240             throw new IllegalArgumentException(message);
 241         }
 242
 243         return this.elevationModel != null ? this.elevationModel.getElevation(latitude, longitude) : 0;
 244     }
 245
 246     public final Vec4 computePointFromPosition(Position position)
 247     {
 248         if (position == null)
 249         {
 250             String message = Logging.getMessage("nullValue.PositionIsNull");
 251             Logging.logger().severe(message);
 252             throw new IllegalArgumentException(message);
 253         }
 254
 255         return this.geodeticToCartesian(position.getLatitude(), position.getLongitude(), position.getElevation());
 256     }
 257
 258     public final Vec4 computePointFromPosition(Angle latitude, Angle longitude, double metersElevation)
 259     {
 260         if (latitude == null || longitude == null)
 261         {
 262             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 263             Logging.logger().severe(message);
 264             throw new IllegalArgumentException(message);
 265         }
 266
 267         return this.geodeticToCartesian(latitude, longitude, metersElevation);
 268     }
 269
 270     public final Position computePositionFromPoint(Vec4 point)
 271     {
 272         if (point == null)
 273         {
 274             String message = Logging.getMessage("nullValue.PointIsNull");
 275             Logging.logger().severe(message);
 276             throw new IllegalArgumentException(message);
 277         }
 278
 279         return this.cartesianToGeodetic(point);
 280     }
 281
 282     public final Position getIntersectionPosition(Line line)
 283     {
 284         if (line == null)
 285         {
 286             String msg = Logging.getMessage("nullValue.LineIsNull");
 287             Logging.logger().severe(msg);
 288             throw new IllegalArgumentException(msg);
 289         }
 290
 291         Intersection[] intersections = this.intersect(line);
 292         if (intersections == null)
 293             return null;
 294
 295         return this.computePositionFromPoint(intersections[0].getIntersectionPoint());
 296     }
 297
 298     // The code below maps latitude / longitude position to globe-centered Cartesian coordinates.
 299     // The Y axis points to the north pole. The Z axis points to the intersection of the prime
 300     // meridian and the equator, in the equatorial plane. The X axis completes a right-handed
 301     // coordinate system, and is 90 degrees east of the Z axis and also in the equatorial plane.
 302
 303     private Vec4 geodeticToCartesian(Angle latitude, Angle longitude, double metersElevation)
 304     {
 305         if (latitude == null || longitude == null)
 306         {
 307             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 308             Logging.logger().severe(message);
 309             throw new IllegalArgumentException(message);
 310         }
 311
 312         double cosLat = latitude.cos();
 313         double sinLat = latitude.sin();
 314
 315         double rpm = // getRadius (in meters) of vertical in prime meridian
 316             this.equatorialRadius / Math.sqrt(1.0 - this.es * sinLat * sinLat);
 317
 318         double x = (rpm + metersElevation) * cosLat * longitude.sin();
 319         double y = (rpm * (1.0 - this.es) + metersElevation) * sinLat;
 320         double z = (rpm + metersElevation) * cosLat * longitude.cos();
 321
 322         return new Vec4(x, y, z);
 323     }
 324
 325     private Position cartesianToGeodetic(Vec4 cart)
 326     {
 327         if (cart == null)
 328         {
 329             String message = Logging.getMessage("nullValue.PointIsNull");
 330             Logging.logger().severe(message);
 331             throw new IllegalArgumentException(message);
 332         }
 333
 334         // according to
 335         // H. Vermeille,
 336         // Direct transformation from geocentric to geodetic ccordinates,
 337         // Journal of Geodesy (2002) 76:451-454
 338         double ra2 = 1 / (this.equatorialRadius * equatorialRadius);
 339
 340         double X = cart.z;
 341         //noinspection SuspiciousNameCombination
 342         double Y = cart.x;
 343         double Z = cart.y;
 344         double e2 = this.es;
 345         double e4 = e2 * e2;
 346
 347         double XXpYY = X * X + Y * Y;
 348         double sqrtXXpYY = Math.sqrt(XXpYY);
 349         double p = XXpYY * ra2;
 350         double q = Z * Z * (1 - e2) * ra2;
 351         double r = 1 / 6.0 * (p + q - e4);
 352         double s = e4 * p * q / (4 * r * r * r);
 353         double t = Math.pow(1 + s + Math.sqrt(s * (2 + s)), 1 / 3.0);
 354         double u = r * (1 + t + 1 / t);
 355         double v = Math.sqrt(u * u + e4 * q);
 356         double w = e2 * (u + v - q) / (2 * v);
 357         double k = Math.sqrt(u + v + w * w) - w;
 358         double D = k * sqrtXXpYY / (k + e2);
 359         double lon = 2 * Math.atan2(Y, X + sqrtXXpYY);
 360         double sqrtDDpZZ = Math.sqrt(D * D + Z * Z);
 361         double lat = 2 * Math.atan2(Z, D + sqrtDDpZZ);
 362         double elevation = (k + e2 - 1) * sqrtDDpZZ / k;
 363
 364         return Position.fromRadians(lat, lon, elevation);
 365     }
 366
 367     public boolean equals(Object o)
 368     {
 369         if (this == o)
 370             return true;
 371         if (o == null || getClass() != o.getClass())
 372             return false;
 373
 374         EllipsoidalGlobe that = (EllipsoidalGlobe) o;
 375
 376         if (Double.compare(that.equatorialRadius, equatorialRadius) != 0)
 377             return false;
 378         if (Double.compare(that.es, es) != 0)
 379             return false;
 380         if (Double.compare(that.polarRadius, polarRadius) != 0)
 381             return false;
 382         if (center != null ? !center.equals(that.center) : that.center != null)
 383             return false;
 384         //noinspection RedundantIfStatement
 385         if (elevationModel != null ? !elevationModel.equals(that.elevationModel) : that.elevationModel != null)
 386             return false;
 387
 388         return true;
 389     }
 390
 391     public int hashCode()
 392     {
 393         int result;
 394         long temp;
 395         temp = equatorialRadius != +0.0d ? Double.doubleToLongBits(equatorialRadius) : 0L;
 396         result = (int) (temp ^ (temp >>> 32));
 397         temp = polarRadius != +0.0d ? Double.doubleToLongBits(polarRadius) : 0L;
 398         result = 31 * result + (int) (temp ^ (temp >>> 32));
 399         temp = es != +0.0d ? Double.doubleToLongBits(es) : 0L;
 400         result = 31 * result + (int) (temp ^ (temp >>> 32));
 401         result = 31 * result + (center != null ? center.hashCode() : 0);
 402         result = 31 * result + (elevationModel != null ? elevationModel.hashCode() : 0);
 403         return result;
 404     }
 405 }