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 2835 2007-09-13 19:29:22Z dcollins $
  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.getMaxElevation();
 104     }
 105
 106     public double getMinElevation()
 107     {
 108         return this.elevationModel.getMinElevation();
 109     }
 110
 111     public double getMaxElevation(Sector sector)
 112     {
 113         return this.elevationModel.getMaxElevation(sector);
 114     }
 115
 116     public double getMinElevation(Sector sector)
 117     {
 118         return this.elevationModel.getMinElevation(sector);
 119     }
 120
 121     public final Extent getExtent()
 122     {
 123         return this;
 124     }
 125
 126     public boolean intersects(Frustum frustum)
 127     {
 128         if (frustum == null)
 129         {
 130             String message = Logging.getMessage("nullValue.FrustumIsNull");
 131             Logging.logger().severe(message);
 132             throw new IllegalArgumentException(message);
 133         }
 134
 135         return frustum.intersects(this);
 136     }
 137
 138     public Intersection[] intersect(Line line)
 139     {
 140         return this.intersect(line, this.equatorialRadius, this.polarRadius);
 141     }
 142
 143     public Intersection[] intersect(Line line, double altitude)
 144     {
 145         return this.intersect(line, this.equatorialRadius + altitude, this.polarRadius + altitude);
 146     }
 147
 148     private Intersection[] intersect(Line line, double equRadius, double polRadius)
 149     {
 150         if (line == null)
 151         {
 152             String message = Logging.getMessage("nullValue.LineIsNull");
 153             Logging.logger().severe(message);
 154             throw new IllegalArgumentException(message);
 155         }
 156
 157         // Taken from Lengyel, 2Ed., Section 5.2.3, page 148.
 158
 159         double m = equRadius / polRadius; // "ratio of the x semi-axis length to the y semi-axis length"
 160         double n = 1d;                    // "ratio of the x semi-axis length to the z semi-axis length"
 161         double m2 = m * m;
 162         double n2 = n * n;
 163         double r2 = equRadius * equRadius; // nominal radius squared //equRadius * polRadius;
 164
 165         double vx = line.getDirection().x;
 166         double vy = line.getDirection().y;
 167         double vz = line.getDirection().z;
 168         double sx = line.getOrigin().x;
 169         double sy = line.getOrigin().y;
 170         double sz = line.getOrigin().z;
 171
 172         double a = vx * vx + m2 * vy * vy + n2 * vz * vz;
 173         double b = 2d * (sx * vx + m2 * sy * vy + n2 * sz * vz);
 174         double c = sx * sx + m2 * sy * sy + n2 * sz * sz - r2;
 175
 176         double discriminant = discriminant(a, b, c);
 177         if (discriminant < 0)
 178             return null;
 179
 180         double discriminantRoot = Math.sqrt(discriminant);
 181         if (discriminant == 0)
 182         {
 183             Vec4 p = line.getPointAt((-b - discriminantRoot) / (2 * a));
 184             return new Intersection[] {new Intersection(p, true)};
 185         }
 186         else // (discriminant > 0)
 187         {
 188             Vec4 near = line.getPointAt((-b - discriminantRoot) / (2 * a));
 189             Vec4 far = line.getPointAt((-b + discriminantRoot) / (2 * a));
 190             return new Intersection[] {new Intersection(near, false), new Intersection(far, false)};
 191         }
 192     }
 193
 194     static private double discriminant(double a, double b, double c)
 195     {
 196         return b * b - 4 * a * c;
 197     }
 198
 199     public boolean intersects(Line line)
 200     {
 201         if (line == null)
 202         {
 203             String msg = Logging.getMessage("nullValue.LineIsNull");
 204             Logging.logger().severe(msg);
 205             throw new IllegalArgumentException(msg);
 206         }
 207
 208         return line.distanceTo(this.center) <= this.equatorialRadius;
 209     }
 210
 211     public boolean intersects(Plane plane)
 212     {
 213         if (plane == null)
 214         {
 215             String msg = Logging.getMessage("nullValue.PlaneIsNull");
 216             Logging.logger().severe(msg);
 217             throw new IllegalArgumentException(msg);
 218         }
 219
 220         double dq1 = plane.dot(this.center);
 221         return dq1 <= this.equatorialRadius;
 222     }
 223
 224     public Vec4 computeSurfaceNormalAtPoint(Vec4 p)
 225     {
 226         if (p == null)
 227         {
 228             String msg = Logging.getMessage("nullValue.PointIsNull");
 229             Logging.logger().severe(msg);
 230             throw new IllegalArgumentException(msg);
 231         }
 232
 233         p = p.subtract3(this.center);
 234
 235         return new Vec4(
 236             p.x / (this.equatorialRadius * this.equatorialRadius),
 237             p.y / (this.polarRadius * this.polarRadius),
 238             p.z / (this.equatorialRadius * this.equatorialRadius)).normalize3();
 239     }
 240
 241     public final ElevationModel getElevationModel()
 242     {
 243         return this.elevationModel;
 244     }
 245
 246     public final double getElevation(Angle latitude, Angle longitude)
 247     {
 248         if (latitude == null || longitude == null)
 249         {
 250             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 251             Logging.logger().severe(message);
 252             throw new IllegalArgumentException(message);
 253         }
 254
 255         return this.elevationModel != null ? this.elevationModel.getElevation(latitude, longitude) : 0;
 256     }
 257
 258     public final Vec4 computePointFromPosition(Position position)
 259     {
 260         if (position == null)
 261         {
 262             String message = Logging.getMessage("nullValue.PositionIsNull");
 263             Logging.logger().severe(message);
 264             throw new IllegalArgumentException(message);
 265         }
 266
 267         return this.geodeticToCartesian(position.getLatitude(), position.getLongitude(), position.getElevation());
 268     }
 269
 270     public final Vec4 computePointFromPosition(Angle latitude, Angle longitude, double metersElevation)
 271     {
 272         if (latitude == null || longitude == null)
 273         {
 274             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 275             Logging.logger().severe(message);
 276             throw new IllegalArgumentException(message);
 277         }
 278
 279         return this.geodeticToCartesian(latitude, longitude, metersElevation);
 280     }
 281
 282     public final Position computePositionFromPoint(Vec4 point)
 283     {
 284         if (point == null)
 285         {
 286             String message = Logging.getMessage("nullValue.PointIsNull");
 287             Logging.logger().severe(message);
 288             throw new IllegalArgumentException(message);
 289         }
 290
 291         return this.cartesianToGeodetic(point);
 292     }
 293
 294     public final Position getIntersectionPosition(Line line)
 295     {
 296         if (line == null)
 297         {
 298             String msg = Logging.getMessage("nullValue.LineIsNull");
 299             Logging.logger().severe(msg);
 300             throw new IllegalArgumentException(msg);
 301         }
 302
 303         Intersection[] intersections = this.intersect(line);
 304         if (intersections == null)
 305             return null;
 306
 307         return this.computePositionFromPoint(intersections[0].getIntersectionPoint());
 308     }
 309
 310     // The code below maps latitude / longitude position to globe-centered Cartesian coordinates.
 311     // The Y axis points to the north pole. The Z axis points to the intersection of the prime
 312     // meridian and the equator, in the equatorial plane. The X axis completes a right-handed
 313     // coordinate system, and is 90 degrees east of the Z axis and also in the equatorial plane.
 314
 315     private Vec4 geodeticToCartesian(Angle latitude, Angle longitude, double metersElevation)
 316     {
 317         if (latitude == null || longitude == null)
 318         {
 319             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 320             Logging.logger().severe(message);
 321             throw new IllegalArgumentException(message);
 322         }
 323
 324         double cosLat = latitude.cos();
 325         double sinLat = latitude.sin();
 326
 327         double rpm = // getRadius (in meters) of vertical in prime meridian
 328             this.equatorialRadius / Math.sqrt(1.0 - this.es * sinLat * sinLat);
 329
 330         double x = (rpm + metersElevation) * cosLat * longitude.sin();
 331         double y = (rpm * (1.0 - this.es) + metersElevation) * sinLat;
 332         double z = (rpm + metersElevation) * cosLat * longitude.cos();
 333
 334         return new Vec4(x, y, z);
 335     }
 336
 337     private Position cartesianToGeodetic(Vec4 cart)
 338     {
 339         if (cart == null)
 340         {
 341             String message = Logging.getMessage("nullValue.PointIsNull");
 342             Logging.logger().severe(message);
 343             throw new IllegalArgumentException(message);
 344         }
 345
 346         // according to
 347         // H. Vermeille,
 348         // Direct transformation from geocentric to geodetic ccordinates,
 349         // Journal of Geodesy (2002) 76:451-454
 350         double ra2 = 1 / (this.equatorialRadius * equatorialRadius);
 351
 352         double X = cart.z;
 353         //noinspection SuspiciousNameCombination
 354         double Y = cart.x;
 355         double Z = cart.y;
 356         double e2 = this.es;
 357         double e4 = e2 * e2;
 358
 359         double XXpYY = X * X + Y * Y;
 360         double sqrtXXpYY = Math.sqrt(XXpYY);
 361         double p = XXpYY * ra2;
 362         double q = Z * Z * (1 - e2) * ra2;
 363         double r = 1 / 6.0 * (p + q - e4);
 364         double s = e4 * p * q / (4 * r * r * r);
 365         double t = Math.pow(1 + s + Math.sqrt(s * (2 + s)), 1 / 3.0);
 366         double u = r * (1 + t + 1 / t);
 367         double v = Math.sqrt(u * u + e4 * q);
 368         double w = e2 * (u + v - q) / (2 * v);
 369         double k = Math.sqrt(u + v + w * w) - w;
 370         double D = k * sqrtXXpYY / (k + e2);
 371         double lon = 2 * Math.atan2(Y, X + sqrtXXpYY);
 372         double sqrtDDpZZ = Math.sqrt(D * D + Z * Z);
 373         double lat = 2 * Math.atan2(Z, D + sqrtDDpZZ);
 374         double elevation = (k + e2 - 1) * sqrtDDpZZ / k;
 375
 376         return Position.fromRadians(lat, lon, elevation);
 377     }
 378
 379     public boolean equals(Object o)
 380     {
 381         if (this == o)
 382             return true;
 383         if (o == null || getClass() != o.getClass())
 384             return false;
 385
 386         EllipsoidalGlobe that = (EllipsoidalGlobe) o;
 387
 388         if (Double.compare(that.equatorialRadius, equatorialRadius) != 0)
 389             return false;
 390         if (Double.compare(that.es, es) != 0)
 391             return false;
 392         if (Double.compare(that.polarRadius, polarRadius) != 0)
 393             return false;
 394         if (center != null ? !center.equals(that.center) : that.center != null)
 395             return false;
 396         //noinspection RedundantIfStatement
 397         if (elevationModel != null ? !elevationModel.equals(that.elevationModel) : that.elevationModel != null)
 398             return false;
 399
 400         return true;
 401     }
 402
 403     public int hashCode()
 404     {
 405         int result;
 406         long temp;
 407         temp = equatorialRadius != +0.0d ? Double.doubleToLongBits(equatorialRadius) : 0L;
 408         result = (int) (temp ^ (temp >>> 32));
 409         temp = polarRadius != +0.0d ? Double.doubleToLongBits(polarRadius) : 0L;
 410         result = 31 * result + (int) (temp ^ (temp >>> 32));
 411         temp = es != +0.0d ? Double.doubleToLongBits(es) : 0L;
 412         result = 31 * result + (int) (temp ^ (temp >>> 32));
 413         result = 31 * result + (center != null ? center.hashCode() : 0);
 414         result = 31 * result + (elevationModel != null ? elevationModel.hashCode() : 0);
 415         return result;
 416     }
 417 }