/*
 * Copyright 2006 Jerry Huxtable
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 * http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package ucar.unidata.geoloc.projection.proj4;

import ucar.unidata.geoloc.ProjectionPoint;
import ucar.unidata.geoloc.ProjectionPointImpl;
import ucar.unidata.geoloc.ProjectionRect;

/**
 * Taken from com.jhlabs.map.proj
 *
 * Also see "http://search.cpan.org/src/DSTAHLKE/Cartography-Projection-GCTP-0.03/gctpc/" for C code
 *
 * @see "http://www.jhlabs.com/java/maps/proj/index.html"
 * @see "http://trac.osgeo.org/proj/"
 *
 * @since Oct 8, 2009
 */

public class MapMath {

  public static final double HALFPI = Math.PI / 2.0;
  public static final double QUARTERPI = Math.PI / 4.0;
  public static final double TWOPI = Math.PI * 2.0;

  public static final ProjectionRect WORLD_BOUNDS_RAD =
      new ProjectionRect(-Math.PI, -Math.PI / 2, Math.PI * 2, Math.PI);
  public static final ProjectionRect WORLD_BOUNDS = new ProjectionRect(-180, -90, 360, 180);

  public static final double EPS10 = 1e-10;
  public static final double RTD = 180.0 / Math.PI;
  public static final double DTR = Math.PI / 180.0;

  /**
   * Degree versions of trigonometric functions
   */
  public static double sind(double v) {
    return Math.sin(v * DTR);
  }

  public static double cosd(double v) {
    return Math.cos(v * DTR);
  }

  public static double tand(double v) {
    return Math.tan(v * DTR);
  }

  public static double asind(double v) {
    return Math.asin(v) * RTD;
  }

  public static double acosd(double v) {
    return Math.acos(v) * RTD;
  }

  public static double atand(double v) {
    return Math.atan(v) * RTD;
  }

  public static double atan2d(double y, double x) {
    return Math.atan2(y, x) * RTD;
  }

  public static double asin(double v) {
    if (Math.abs(v) > 1.)
      return v < 0.0 ? -Math.PI / 2 : Math.PI / 2;
    return Math.asin(v);
  }

  public static double acos(double v) {
    if (Math.abs(v) > 1.)
      return v < 0.0 ? Math.PI : 0.0;
    return Math.acos(v);
  }

  public static double sqrt(double v) {
    return v < 0.0 ? 0.0 : Math.sqrt(v);
  }

  public static double distance(double dx, double dy) {
    return Math.sqrt(dx * dx + dy * dy);
  }

  public static double distance(ProjectionPoint a, ProjectionPoint b) {
    return distance(a.getX() - b.getX(), a.getY() - b.getY());
  }

  public static double hypot(double x, double y) {
    if (x < 0.0)
      x = -x;
    else if (x == 0.0)
      return y < 0.0 ? -y : y;
    if (y < 0.0)
      y = -y;
    else if (y == 0.0)
      return x;
    if (x < y) {
      x /= y;
      return y * Math.sqrt(1.0 + x * x);
    } else {
      y /= x;
      return x * Math.sqrt(1.0 + y * y);
    }
  }

  public static double atan2(double y, double x) {
    return Math.atan2(y, x);
  }

  public static double trunc(double v) {
    return v < 0.0 ? Math.ceil(v) : Math.floor(v);
  }

  public static double frac(double v) {
    return v - trunc(v);
  }

  public static double degToRad(double v) {
    return v * Math.PI / 180.0;
  }

  public static double radToDeg(double v) {
    return v * 180.0 / Math.PI;
  }

  // For negative angles, d should be negative, m & s positive.
  public static double dmsToRad(double d, double m, double s) {
    if (d >= 0)
      return (d + m / 60 + s / 3600) * Math.PI / 180.0;
    return (d - m / 60 - s / 3600) * Math.PI / 180.0;
  }

  // For negative angles, d should be negative, m & s positive.
  public static double dmsToDeg(double d, double m, double s) {
    if (d >= 0)
      return (d + m / 60 + s / 3600);
    return (d - m / 60 - s / 3600);
  }

  public static double normalizeLatitude(double angle) {
    if (Double.isInfinite(angle) || Double.isNaN(angle))
      throw new RuntimeException("Infinite latitude");

    while (angle > MapMath.HALFPI)
      angle -= Math.PI;
    while (angle < -MapMath.HALFPI)
      angle += Math.PI;
    return angle;
    // return Math.IEEEremainder(angle, Math.PI);
  }

  public static double normalizeLongitude(double angle) {
    if (Double.isInfinite(angle) || Double.isNaN(angle))
      throw new RuntimeException("Infinite longitude");

    while (angle > Math.PI)
      angle -= TWOPI;
    while (angle < -Math.PI)
      angle += TWOPI;
    return angle;
    // return Math.IEEEremainder(angle, Math.PI);
  }

  public static double normalizeAngle(double angle) {
    if (Double.isInfinite(angle) || Double.isNaN(angle))
      throw new RuntimeException("Infinite angle");

    while (angle > TWOPI)
      angle -= TWOPI;
    while (angle < 0)
      angle += TWOPI;
    return angle;
  }

  /*
   * public static void latLongToXYZ(ProjectionPoint ll, Point3D xyz) {
   * double c = Math.cos(ll.y);
   * xyz.x = c * Math.cos(ll.x);
   * xyz.y = c * Math.sin(ll.x);
   * xyz.z = Math.sin(ll.y);
   * }
   * 
   * public static void xyzToLatLong(Point3D xyz, ProjectionPoint ll) {
   * ll.y = MapMath.asin(xyz.z);
   * ll.x = MapMath.atan2(xyz.y, xyz.x);
   * }
   */

  public static double greatCircleDistance(double lon1, double lat1, double lon2, double lat2) {
    double dlat = Math.sin((lat2 - lat1) / 2);
    double dlon = Math.sin((lon2 - lon1) / 2);
    double r = Math.sqrt(dlat * dlat + Math.cos(lat1) * Math.cos(lat2) * dlon * dlon);
    return 2.0 * Math.asin(r);
  }

  public static double sphericalAzimuth(double lat0, double lon0, double lat, double lon) {
    double diff = lon - lon0;
    double coslat = Math.cos(lat);

    return Math.atan2(coslat * Math.sin(diff),
        (Math.cos(lat0) * Math.sin(lat) - Math.sin(lat0) * coslat * Math.cos(diff)));
  }

  public static boolean sameSigns(double a, double b) {
    return a < 0 == b < 0;
  }

  public static boolean sameSigns(int a, int b) {
    return a < 0 == b < 0;
  }

  public static double takeSign(double a, double b) {
    a = Math.abs(a);
    if (b < 0)
      return -a;
    return a;
  }

  public static int takeSign(int a, int b) {
    a = Math.abs(a);
    if (b < 0)
      return -a;
    return a;
  }

  public static final int DONT_INTERSECT = 0;
  public static final int DO_INTERSECT = 1;
  public static final int COLLINEAR = 2;

  public static int intersectSegments(ProjectionPoint aStart, ProjectionPoint aEnd, ProjectionPoint bStart,
      ProjectionPoint bEnd, ProjectionPointImpl p) {
    double a1, a2, b1, b2, c1, c2;
    double r1, r2, r3, r4;
    double denom, offset, num;

    a1 = aEnd.getY() - aStart.getY();
    b1 = aStart.getX() - aEnd.getX();
    c1 = aEnd.getX() * aStart.getY() - aStart.getX() * aEnd.getY();
    r3 = a1 * bStart.getX() + b1 * bStart.getY() + c1;
    r4 = a1 * bEnd.getX() + b1 * bEnd.getY() + c1;

    if (r3 != 0 && r4 != 0 && sameSigns(r3, r4))
      return DONT_INTERSECT;

    a2 = bEnd.getY() - bStart.getY();
    b2 = bStart.getX() - bEnd.getX();
    c2 = bEnd.getX() * bStart.getY() - bStart.getX() * bEnd.getY();
    r1 = a2 * aStart.getX() + b2 * aStart.getY() + c2;
    r2 = a2 * aEnd.getX() + b2 * aEnd.getY() + c2;

    if (r1 != 0 && r2 != 0 && sameSigns(r1, r2))
      return DONT_INTERSECT;

    denom = a1 * b2 - a2 * b1;
    if (denom == 0)
      return COLLINEAR;

    offset = denom < 0 ? -denom / 2 : denom / 2;

    num = b1 * c2 - b2 * c1;
    p.setX((num < 0 ? num - offset : num + offset) / denom);

    num = a2 * c1 - a1 * c2;
    p.setY((num < 0 ? num - offset : num + offset) / denom);

    return DO_INTERSECT;
  }

  public static double dot(ProjectionPoint a, ProjectionPoint b) {
    return a.getX() * b.getX() + a.getY() * b.getY();
  }

  public static ProjectionPoint perpendicular(ProjectionPoint a) {
    return ProjectionPoint.create(-a.getY(), a.getX());
  }

  public static ProjectionPoint add(ProjectionPoint a, ProjectionPoint b) {
    return ProjectionPoint.create(a.getX() + b.getX(), a.getY() + b.getY());
  }

  public static ProjectionPoint subtract(ProjectionPoint a, ProjectionPoint b) {
    return ProjectionPoint.create(a.getX() - b.getX(), a.getY() - b.getY());
  }

  public static ProjectionPoint multiply(ProjectionPoint a, ProjectionPoint b) {
    return ProjectionPoint.create(a.getX() * b.getX(), a.getY() * b.getY());
  }

  public static double cross(ProjectionPoint a, ProjectionPoint b) {
    return a.getX() * b.getY() - b.getX() * a.getY();
  }

  public static double cross(double x1, double y1, double x2, double y2) {
    return x1 * y2 - x2 * y1;
  }

  public static void normalize(ProjectionPointImpl a) {
    double d = distance(a.getX(), a.getY());
    a.setLocation(a.getX() / d, a.getY() / d);
  }

  public static void negate(ProjectionPointImpl a) {
    a.setLocation(-a.getX(), -a.getY());
  }

  public static double longitudeDistance(double l1, double l2) {
    return Math.min(Math.abs(l1 - l2),
        ((l1 < 0) ? l1 + Math.PI : Math.PI - l1) + ((l2 < 0) ? l2 + Math.PI : Math.PI - l2));
  }

  public static double geocentricLatitude(double lat, double flatness) {
    double f = 1.0 - flatness;
    return Math.atan((f * f) * Math.tan(lat));
  }

  public static double geographicLatitude(double lat, double flatness) {
    double f = 1.0 - flatness;
    return Math.atan(Math.tan(lat) / (f * f));
  }

  public static double tsfn(double phi, double sinphi, double e) {
    sinphi *= e;
    return (Math.tan(.5 * (MapMath.HALFPI - phi)) / Math.pow((1. - sinphi) / (1. + sinphi), .5 * e));
  }

  public static double msfn(double sinphi, double cosphi, double es) {
    return cosphi / Math.sqrt(1.0 - es * sinphi * sinphi);
  }

  private static final int N_ITER = 15;

  public static double phi2(double ts, double e) {
    double eccnth, phi, con, dphi;
    int i;

    eccnth = .5 * e;
    phi = MapMath.HALFPI - 2. * Math.atan(ts);
    i = N_ITER;
    do {
      con = e * Math.sin(phi);
      dphi = MapMath.HALFPI - 2. * Math.atan(ts * Math.pow((1. - con) / (1. + con), eccnth)) - phi;
      phi += dphi;
    } while (Math.abs(dphi) > 1e-10 && --i != 0);
    if (i <= 0)
      throw new RuntimeException();
    return phi;
  }

  private static final double C00 = 1.0;
  private static final double C02 = .25;
  private static final double C04 = .046875;
  private static final double C06 = .01953125;
  private static final double C08 = .01068115234375;
  private static final double C22 = .75;
  private static final double C44 = .46875;
  private static final double C46 = .01302083333333333333;
  private static final double C48 = .00712076822916666666;
  private static final double C66 = .36458333333333333333;
  private static final double C68 = .00569661458333333333;
  private static final double C88 = .3076171875;
  private static final int MAX_ITER = 10;

  public static double[] enfn(double es) {
    double t;
    double[] en = new double[5];
    en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
    en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
    en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
    en[3] = (t *= es) * (C66 - es * C68);
    en[4] = t * es * C88;
    return en;
  }

  public static double mlfn(double phi, double sphi, double cphi, double[] en) {
    cphi *= sphi;
    sphi *= sphi;
    return en[0] * phi - cphi * (en[1] + sphi * (en[2] + sphi * (en[3] + sphi * en[4])));
  }

  public static double inv_mlfn(double arg, double es, double[] en) {
    double s, t, phi, k = 1. / (1. - es);

    phi = arg;
    for (int i = MAX_ITER; i != 0; i--) {
      s = Math.sin(phi);
      t = 1. - es * s * s;
      phi -= t = (mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
      if (Math.abs(t) < 1e-11)
        return phi;
    }
    return phi;
  }

  private static final double P00 = .33333333333333333333;
  private static final double P01 = .17222222222222222222;
  private static final double P02 = .10257936507936507936;
  private static final double P10 = .06388888888888888888;
  private static final double P11 = .06640211640211640211;
  private static final double P20 = .01641501294219154443;

  public static double[] authset(double es) {
    double t;
    double[] APA = new double[3];
    APA[0] = es * P00;
    t = es * es;
    APA[0] += t * P01;
    APA[1] = t * P10;
    t *= es;
    APA[0] += t * P02;
    APA[1] += t * P11;
    APA[2] = t * P20;
    return APA;
  }

  public static double authlat(double beta, double[] APA) {
    double t = beta + beta;
    return (beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t + t) + APA[2] * Math.sin(t + t + t));
  }

  public static double qsfn(double sinphi, double e, double one_es) {
    double con;

    if (e >= 1.0e-7) {
      con = e * sinphi;
      return (one_es * (sinphi / (1. - con * con) - (.5 / e) * Math.log((1. - con) / (1. + con))));
    } else
      return (sinphi + sinphi);
  }

  /*
   * Java translation of "Nice Numbers for Graph Labels"
   * by Paul Heckbert
   * from "Graphics Gems", Academic Press, 1990
   */
  public static double niceNumber(double x, boolean round) {
    int expv; /* exponent of x */
    double f; /* fractional part of x */
    double nf; /* nice, rounded fraction */

    expv = (int) Math.floor(Math.log(x) / Math.log(10));
    f = x / Math.pow(10., expv); /* between 1 and 10 */
    if (round) {
      if (f < 1.5)
        nf = 1.;
      else if (f < 3.)
        nf = 2.;
      else if (f < 7.)
        nf = 5.;
      else
        nf = 10.;
    } else if (f <= 1.)
      nf = 1.;
    else if (f <= 2.)
      nf = 2.;
    else if (f <= 5.)
      nf = 5.;
    else
      nf = 10.;
    return nf * Math.pow(10., expv);
  }
}
