/*
Copyright (C) 2001, 2006 United States Government
as represented by the Administrator of the
National Aeronautics and Space Administration.
All Rights Reserved.
*/
package gov.nasa.worldwind.globes;

import gov.nasa.worldwind.*;
import gov.nasa.worldwind.geom.*;

/**
 * @author Tom Gaskins
 * @version $Id: EllipsoidalGlobe.java 1782 2007-05-08 06:27:54Z tgaskins $
 */
public class EllipsoidalGlobe extends WWObjectImpl implements Globe
{
    private final double equatorialRadius;
    private final double polarRadius;
    private final double es;
    private final Point center;

    private final ElevationModel elevationModel;

    public EllipsoidalGlobe(double equatorialRadius, double polarRadius, double es, ElevationModel em)
    {
        if (em == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.ElevationModelIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        this.equatorialRadius = equatorialRadius;
        this.polarRadius = polarRadius;
        this.es = es; // assume it's consistent with the two radii
        this.center = Point.ZERO;
        this.elevationModel = em;
    }

    public final double getRadius()
    {
        return this.equatorialRadius;
    }

    public final double getEquatorialRadius()
    {
        return this.equatorialRadius;
    }

    public final double getPolarRadius()
    {
        return this.polarRadius;
    }

    public double getMaximumRadius()
    {
        return this.equatorialRadius;
    }

    public double getRadiusAt(Angle latitude, Angle longitude)
    {
        if (latitude == null || longitude == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.AngleIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        return this.computePointFromPosition(latitude, longitude, 0d).length();
    }

    public double getEccentricitySquared()
    {
        return this.es;
    }

    public final double getDiameter()
    {
        return this.equatorialRadius * 2;
    }

    public final Point getCenter()
    {
        return this.center;
    }

    public double getMaxElevation()
    {
        return this.elevationModel.getMaximumElevation();
    }

    public double getMinElevation()
    {
        return this.elevationModel.getMinimumElevation();
    }

    public final Extent getExtent()
    {
        return this;
    }

    public boolean intersects(Frustum frustum)
    {
        if (frustum == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.FrustumIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        return frustum.intersects(this);
    }

    public Intersection[] intersect(Line line)
    {
        if (line == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.LineIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        // Taken from Lengyel, 2Ed., Section 5.2.3, page 148.
        double m = this.equatorialRadius / this.polarRadius;
        double n = 1d; //this.equatorialRadius / this.equatorialRadius;
        double m2 = m * m;
        double n2 = n * n;

        double vx = line.getDirection().getX();
        double vy = line.getDirection().getY();
        double vz = line.getDirection().getZ();
        double sx = line.getOrigin().getX();
        double sy = line.getOrigin().getY();
        double sz = line.getOrigin().getZ();

        double a = vx * vx + m2 * vy * vy + n2 * vz * vz;
        double b = 2d * (sx * vx + m2 * sy * vy + n2 * sz * vz);
        double c = sx * sx + m2 * sy * sy + n2 * sz * sz - this.equatorialRadius * this.equatorialRadius;

        double discriminant = discriminant(a, b, c);
        if (discriminant < 0)
            return null;

        double discriminantRoot = Math.sqrt(discriminant);
        if (discriminant == 0)
        {
            Point p = line.getPointAt((-b - discriminantRoot) / (2 * a));
            return new Intersection[] {new Intersection(p, true)};
        }
        else // (discriminant > 0)
        {
            Point near = line.getPointAt((-b - discriminantRoot) / (2 * a));
            Point far = line.getPointAt((-b + discriminantRoot) / (2 * a));
            return new Intersection[] {new Intersection(near, false), new Intersection(far, false)};
        }
    }

    static private double discriminant(double a, double b, double c)
    {
        return b * b - 4 * a * c;
    }

    public boolean intersects(Line line)
    {
        if (line == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.LineIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        return line.distanceTo(this.center) <= this.equatorialRadius;
    }

    public boolean intersects(Plane plane)
    {
        if (plane == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.PlaneIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        double dq1 = plane.dot(this.center);
        return dq1 <= this.equatorialRadius;
    }

    public Point computeSurfaceNormalAtPoint(Point p)
    {
        if (p == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.PointIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        p = p.subtract(this.center);

        return new Point(p.x() / (this.equatorialRadius * this.equatorialRadius),
            p.y() / (this.polarRadius * this.polarRadius),
            p.z() / (this.equatorialRadius * this.equatorialRadius)).normalize();
    }

    public final ElevationModel getElevationModel()
    {
        return this.elevationModel;
    }

    public final double getElevation(Angle latitude, Angle longitude)
    {
        if (latitude == null || longitude == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.LatitudeOrLongitudeIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        return this.elevationModel != null ? this.elevationModel.getElevation(latitude, longitude) : 0;
    }

    public final Point computePointFromPosition(Position position)
    {
        if (position == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.PositionIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        return this.geodeticToCartesian(position.getLatitude(), position.getLongitude(), position.getElevation());
    }

    public final Point computePointFromPosition(Angle latitude, Angle longitude, double metersElevation)
    {
        if (latitude == null || longitude == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.LatitudeOrLongitudeIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        return this.geodeticToCartesian(latitude, longitude, metersElevation);
    }

    public final Position computePositionFromPoint(Point point)
    {
        if (point == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.PointIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        return this.cartesianToGeodetic(point);
    }

    public final Position getIntersectionPosition(Line line)
    {
        if (line == null)
        {
            String msg = WorldWind.retrieveErrMsg("nullValue.LineIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, msg);
            throw new IllegalArgumentException(msg);
        }

        Intersection[] intersections = this.intersect(line);
        if (intersections == null)
            return null;

        return this.computePositionFromPoint(intersections[0].getIntersectionPoint());
    }

    // The code below maps latitude / longitude position to globe-centered Cartesian coordinates.
    // The Y axis points to the north pole. The Z axis points to the intersection of the prime
    // meridian and the equator, in the equatorial plane. The X axis completes a right-handed
    // coordinate system, and is 90 degrees east of the Z axis and also in the equatorial plane.

    private Point geodeticToCartesian(Angle latitude, Angle longitude, double metersElevation)
    {
        if (latitude == null || longitude == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.LatitudeOrLongitudeIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        double cosLat = latitude.cos();
        double sinLat = latitude.sin();

        double rpm = // getRadius (in meters) of vertical in prime meridian
            this.equatorialRadius / Math.sqrt(1.0 - this.es * sinLat * sinLat);

        double x = (rpm + metersElevation) * cosLat * longitude.sin();
        double y = (rpm * (1.0 - this.es) + metersElevation) * sinLat;
        double z = (rpm + metersElevation) * cosLat * longitude.cos();

        return new Point(x, y, z);
    }

    private Position cartesianToGeodetic(Point cart)
    {
        if (cart == null)
        {
            String message = WorldWind.retrieveErrMsg("nullValue.PointIsNull");
            WorldWind.logger().log(java.util.logging.Level.FINE, message);
            throw new IllegalArgumentException(message);
        }

        // according to
        // H. Vermeille,
        // Direct transformation from geocentric to geodetic ccordinates,
        // Journal of Geodesy (2002) 76:451-454
        double ra2 = 1 / (this.equatorialRadius * equatorialRadius);

        double X = cart.z();
        double Y = cart.x();
        double Z = cart.y();
        double e2 = this.es;
        double e4 = e2 * e2;

        double XXpYY = X * X + Y * Y;
        double sqrtXXpYY = Math.sqrt(XXpYY);
        double p = XXpYY * ra2;
        double q = Z * Z * (1 - e2) * ra2;
        double r = 1 / 6.0 * (p + q - e4);
        double s = e4 * p * q / (4 * r * r * r);
        double t = Math.pow(1 + s + Math.sqrt(s * (2 + s)), 1 / 3.0);
        double u = r * (1 + t + 1 / t);
        double v = Math.sqrt(u * u + e4 * q);
        double w = e2 * (u + v - q) / (2 * v);
        double k = Math.sqrt(u + v + w * w) - w;
        double D = k * sqrtXXpYY / (k + e2);
        double lon = 2 * Math.atan2(Y, X + sqrtXXpYY);
        double sqrtDDpZZ = Math.sqrt(D * D + Z * Z);
        double lat = 2 * Math.atan2(Z, D + sqrtDDpZZ);
        double elevation = (k + e2 - 1) * sqrtDDpZZ / k;

        return Position.fromRadians(lat, lon, elevation);
    }

    public boolean equals(Object o)
    {
        if (this == o)
            return true;
        if (o == null || getClass() != o.getClass())
            return false;

        EllipsoidalGlobe that = (EllipsoidalGlobe) o;

        if (Double.compare(that.equatorialRadius, equatorialRadius) != 0)
            return false;
        if (Double.compare(that.es, es) != 0)
            return false;
        if (Double.compare(that.polarRadius, polarRadius) != 0)
            return false;
        if (center != null ? !center.equals(that.center) : that.center != null)
            return false;
        //noinspection RedundantIfStatement
        if (elevationModel != null ? !elevationModel.equals(that.elevationModel) : that.elevationModel != null)
            return false;

        return true;
    }

    public int hashCode()
    {
        int result;
        long temp;
        temp = equatorialRadius != +0.0d ? Double.doubleToLongBits(equatorialRadius) : 0L;
        result = (int) (temp ^ (temp >>> 32));
        temp = polarRadius != +0.0d ? Double.doubleToLongBits(polarRadius) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        temp = es != +0.0d ? Double.doubleToLongBits(es) : 0L;
        result = 31 * result + (int) (temp ^ (temp >>> 32));
        result = 31 * result + (center != null ? center.hashCode() : 0);
        result = 31 * result + (elevationModel != null ? elevationModel.hashCode() : 0);
        return result;
    }
}