#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""2D geometry objects.

:Date: 2019-08-28

.. module:: geometry
  :platform: *nix, Windows
  :synopsis: Geometry objects.

.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
"""
import math
import numpy as np


def distance(point1, point2):
  """Distance between two points (or length of a straight line).

  :param point1: first point (first end point of straight line)
  :type point1: tuple
  :param point2: second point (second end point of straight line)
  :type point2: tuple

  :returns: distance between the two points
  :rtype: float
  """
  return math.sqrt((point2[0]-point1[0])**2 + (point2[1]-point1[1])**2)


def angle(point1, point2=None):
  """Angle of point or between two points.

  :param point1: (first) point
  :type point1: tuple
  :param point2: second point (default = None)
  :type point2: tuple

  :returns: angle of point or between two points
  :rtype: float
  """
  if point2 is None:
    return math.atan2(point1[1], point1[0])
  return math.atan2(point2[1]-point1[1], point2[0]-point1[0])


def translate(vec, *pts):
  """Translate a point or polygon by a given vector.

  :param vec: translation vector
  :type vec: tuple
  :param `*pts`: points to translate

  :returns: (point_x, point_y) or (point1, point2, ...)
  :rtype: tuple

  .. seealso::
    :meth:`translate_xy`
  """
  vx, vy = vec
  return tuple([(x+vx, y+vy) for (x, y) in pts])


def translate_xy(vec, x, y):
  """Translate a point or polygon by a given vector.

  :param vec: translation vector
  :type vec: tuple
  :param x: points to translate
  :type x: int or float or list
  :param y: points to translate
  :type y: int or float or list

  :returns: (x', y')
  :rtype: tuple

  .. seealso::
    :meth:`translate`
  """
  vx, vy = vec

  if not hasattr(x, "__len__"):
    x = [x]

  if not hasattr(y, "__len__"):
    y = [y]

  xp = [xi+vx for xi in x]
  yp = [yi+vy for yi in y]

  # no list if it is only one value
  if len(xp) == 1:
    return xp[0], yp[0]

  return xp, yp


def rotate(origin, angle, *pts, **kwargs):
  """Rotate a point or polygon counterclockwise by a given angle
  around a given origin. The angle should be given in radians.

  :param origin: the center of rotation
  :type origin: tuple
  :param angle: the rotation angle
  :type angle: int or float
  :param `*pts`: points to rotate
  :param `**kwargs`: options

  :returns: (point_x, point_y) or (point1, point2, ...)
  :rtype: tuple

  .. seealso::
    :meth:`rotate_xy`
  """
  ox, oy = origin

  # add first point to the end
  if "closed" in kwargs and kwargs["closed"] is True:
    pts += (pts[0],)

  result = tuple([(ox + math.cos(angle) * (px - ox) - math.sin(angle) * (py - oy),
                   oy + math.sin(angle) * (px - ox) + math.cos(angle) * (py - oy))
                  for (px, py) in pts])

  # no tuple in tuple if it is only one point
  if len(pts) == 1:
    return result[0][0], result[0][1]

  return result


def rotate_deg(origin, angle, *pts, **kwargs):
  """Rotate a point or polygon counterclockwise by a given angle
  around a given origin. The angle should be given in degrees.

  :param origin: the center of rotation
  :type origin: tuple
  :param angle: the rotation angle
  :type angle: int or float
  :param `*pts`: points to rotate
  :param `**kwargs`: options

  :returns: (point_x, point_y) or (point1, point2, ...)
  :rtype: tuple

  .. seealso::
    :meth:`rotate`
  """
  return rotate(origin, angle*math.pi/180, *pts, **kwargs)


def rotate_xy(origin, angle, x, y, **kwargs):
  """Rotate x and y coordinates counterclockwise by a given angle
  around a given origin. The angle should be given in radians.

  :param origin: the center of rotation
  :type origin: tuple
  :param angle: the rotation angle
  :type angle: int or float
  :param x: x coordinates
  :type x: int or float or list
  :param y: y coordinates
  :type y: int or float or list
  :param `**kwargs`: options

  .. seealso::
    :meth:`rotate`
  """
  ox, oy = origin

  if not hasattr(x, "__len__"):
    x = [x]

  if not hasattr(y, "__len__"):
    y = [y]

  # add first point to the end
  if "closed" in kwargs and kwargs["closed"] is True:
    x.append(x[0])
    y.append(y[0])

  x_result = [ox + math.cos(angle) * (xi - ox) - math.sin(angle) * (yi - oy)
              for xi, yi in zip(x, y)]

  y_result = [oy + math.sin(angle) * (xi - ox) + math.cos(angle) * (yi - oy)
              for xi, yi in zip(x, y)]

  # no list if it is only one value
  if len(x_result) == 1:
    return x_result[0], y_result[0]

  return x_result, y_result


def rectangle(width, height):
  """\
  :param width: the width of the rectangle
  :type width: int or float
  :param height: the height of the rectangle
  :type height: int or float

  :returns: (point1, point2, point3, point4)
  :rtype: tuple
  """
  pt1 = (-width/2, -height/2)
  pt2 = (width/2, -height/2)
  pt3 = (width/2, height/2)
  pt4 = (-width/2, height/2)
  return pt1, pt2, pt3, pt4, pt1


def square(width):
  """\
  :param width: the edge size of the square
  :type width: int or float

  :returns: (point1, point2, point3, point4)
  :rtype: tuple

  .. seealso::
    :meth:`rectangle`
  """
  return rectangle(width, width)


def lines(pts, **kwargs):
  """Lines defined by a list of end points.

  :param pts: list of points in absolute global coordinate system. If
    keyword inc is given than the inc decides what the left and the
    right end point of the line is, otherwise it is assumed that the
    points build a solid line, that is lines between the given points
    in given order.
  :type pts: list

  :param `**kwargs`: options:

    * deformation -- list of points. Additional deformation
      (translation) at point.
    * factor -- factor of the deformation (default = 1).
    * inc -- the incidence table, a list of 2 element lists. The inc
      decides what the left and the right end point of the line is.
    * index_offset -- starting index of lists (default = 0).

  :returns: list of endpoints for each line;
    [((point1_x, point1_y), (point2_x, point2_y)),
    (p1, p2),
    ...]
  :rtype: list

  .. seealso::
    :meth:`~geometry_plot.plot_lines` of the :mod:`geometry_plot`
    module to plot the lines
  """
  if 'index_offset' not in kwargs:
    kwargs['index_offset'] = 0
  if 'deformation' in kwargs:
    if 'factor' not in kwargs:
      kwargs['factor'] = 1
    pts = [(p[0]+d[0]*kwargs['factor'], p[1]+d[1]*kwargs['factor']) for
           p, d in zip(pts, kwargs['deformation'])]
  if 'inc' in kwargs:
    return [(pts[l-kwargs['index_offset']],
             pts[r-kwargs['index_offset']]) for l, r in kwargs['inc']]
  return list(zip(pts[:-1], pts[1:]))


def cubics(pts, **kwargs):
  """Cubic lines defined by a list of two end points. The deformation
  as displacement and rotation (radians) is defined element wise as
  keyword argument deformation or global node wise as
  global_deformation. The global coordinate system is xy. x in the
  right direction and y in the top direction.

  :param pts: list of points in absolute global coordinate system. If
    keyword inc is given than the inc decides what the left and the
    right end point of the line is, otherwise it is assumed that the
    points build a solid line, that is lines between the given points
    in given order.
  :type pts_rot: list

  :param `**kwargs`: options:

    * deformation -- list of deformation element wise. Additional
      deformation (translation and rotation in radians) at element
      left and right node.
    * rotation_plane -- rotation plane of the element wise
      deformation defined by a string; either 'xy' or 'xz' (default
      = 'xy'). x in the right direction and y in the top direction
      or z in the bottom direction.
    * global_deformation -- list of deformation global node wise.
      Additional deformation (horizontal translation, vertical
      translation and rotation in radians) at node.
    * factor -- factor of the derformation (default = 1).
    * inc -- the incidence table, a list of 2 element lists. The inc
      decides what the left and the right end point of the line is.
    * index_offset -- starting index of lists (default = 0).

  :returns: list of endpoints for each line;
    [(((point1_x, point1_y) angle1), ((point2_x, point2_y), angle2),
    (p1, angle1, p2, angle2),
    ...]
  :rtype: list
  """
  if 'index_offset' not in kwargs:
    kwargs['index_offset'] = 0
  if 'deformation' in kwargs or 'global_deformation' in kwargs:
    if 'factor' not in kwargs:
      kwargs['factor'] = 1
    if 'inc' in kwargs:
      if 'global_deformation' in kwargs:
        lr = [(pts[l-kwargs['index_offset']],
               pts[r-kwargs['index_offset']]) for l, r in kwargs['inc']]
        ang = [angle(l, r) for l, r in lr]
        # system deformation
        U = kwargs['global_deformation']
        if 'rotation_plane' in kwargs and 'xz' == kwargs['rotation_plane']:
          # system deformation left X Z RY right X Z RY element wise
          # Z downwards -> convert to X Y RZ
          Ue = [[ U[(l-kwargs['index_offset'])][0],
                 -U[(l-kwargs['index_offset'])][1],
                  U[(l-kwargs['index_offset'])][2],
                  U[(r-kwargs['index_offset'])][0],
                 -U[(r-kwargs['index_offset'])][1],
                  U[(r-kwargs['index_offset'])][2]] for l, r in kwargs['inc']]
        else:
          # system deformation left X Y RZ right X Y RZ element wise
          Ue = [[U[(l-kwargs['index_offset'])][0],
                 U[(l-kwargs['index_offset'])][1],
                 U[(l-kwargs['index_offset'])][2],
                 U[(r-kwargs['index_offset'])][0],
                 U[(r-kwargs['index_offset'])][1],
                 U[(r-kwargs['index_offset'])][2]] for l, r in kwargs['inc']]
        # element deformation, X Y RZ to x r rz
        # back transformation T^T = [[c, s, 0], [-s, c, 0], [0, 0, 1]]
        u = [[( math.cos(angi)*Uei[0]+math.sin(angi)*Uei[1])*kwargs['factor'],
              (-math.sin(angi)*Uei[0]+math.cos(angi)*Uei[1])*kwargs['factor'],
                               Uei[2]                       *kwargs['factor'],
              ( math.cos(angi)*Uei[3]+math.sin(angi)*Uei[4])*kwargs['factor'],
              (-math.sin(angi)*Uei[3]+math.cos(angi)*Uei[4])*kwargs['factor'],
                               Uei[5]                       *kwargs['factor']]
        for Uei, angi in zip(Ue, ang)]
      else:  # deformation
        # the deformation is in element coordinate system, therefore the angle is needed
        if 'rotation_plane' in kwargs and 'xz' == kwargs['rotation_plane']:
          u = [[ ue[0]*kwargs['factor'],
                -ue[1]*kwargs['factor'],
                 ue[2]*kwargs['factor'],
                 ue[3]*kwargs['factor'],
                -ue[4]*kwargs['factor'],
                 ue[5]*kwargs['factor']] for ue in kwargs['deformation']]
        else:
          u = [[ui*kwargs['factor'] for ui in ue] for ue in kwargs['deformation']]
      return [(pts[l-kwargs['index_offset']],
               pts[r-kwargs['index_offset']],
               d) for (l, r), d in zip(kwargs['inc'], u)]
  return list(zip(pts[:-1], pts[1:]))


def interpolate_hermite(lvd, lr, rvd, rr, lhd=0, rhd=0, scale_x=1, scale_y=1, samples=10):
  r"""Interpolate cubic line with hermite boundary conditions.

  :param lvd: left vertcal deflection
  :type lvd: int or float
  :param lr: left rotation
  :type lr: int or float
  :param rvd: right vertical deflection
  :type rvd: int or float
  :param rr: right rotation
  :type rr: int or float
  :param lhd: left horizontal deformation (default = 0)
  :type lhd: int or float
  :param rhd: right horizontal deformation (default = 0)
  :type rhd: int or float
  :param scale_x: length of element (default = 1)
  :type scale_x: int or float
  :param scale_y: factor of the deformation (default = 1).
    This does not change the length.
  :type scale_y: int or float
  :param samples: number of sampling points (default = 10)
  :type samples: int

  .. math::
    s = \frac{x - x_1}{L} \\
    x = s\,L + x_1

  """
  L = scale_x
  Lp = L + rhd - lhd

  # x=[0,1] in non-dimensional coordinates
  x = np.linspace(0, 1, num=samples)
  N1 =  1     - 3*x**2 + 2*x**3
  N2 = (    x - 2*x**2 +   x**3)*Lp
  N3 =          3*x**2 - 2*x**3
  N4 = (      -   x**2 +   x**3)*Lp
  x = x*Lp


  # x=[0,L] in global coordinates
  x = np.linspace(0, Lp, num=samples)
  N1 = 1     - 3*x**2/Lp**2 + 2*x**3/Lp**3
  N2 =     x - 2*x**2/Lp    +   x**3/Lp**2
  N3 =         3*x**2/Lp**2 - 2*x**3/Lp**3
  N4 =       -   x**2/Lp    +   x**3/Lp**2


  v = N1*lvd + N2*lr + N3*rvd + N4*rr
  x = x + lhd
  y = v*scale_y

  return x, y


#
# matplotlib format, return lists for x and y
#

def line(point1, point2, samples=2):
  """Line defined by two end points.

  .. math::
    y = \\frac{y_2-y_1}{x_2-x_1}(x-x_1) + y_1

  :param point1: one end point
  :type point1: tuple
  :param point2: other end point
  :type point2: tuple
  :param samples: number of sampling points (default = 2)
  :type samples: int

  :returns: ((point1_x, point2_x), (points1_y, point2_y)) or
    ([sample_point1_x, sample_point2_x, ...],
    [sample_points1_y, sample_point2_y, ...])
  :rtype: tuple

  :Example:

  >>> x, y = line((0, 0), (1, 0))
  >>> print(x, y)
  ((0, 1), (0, 0))
  """
  p1x, p1y = point1
  p2x, p2y = point2

  denominator = (p1x - p2x)

  if samples > 2 and denominator > 0:
    x = np.linspace(p1x, p2x, samples)
    a = (p1y - p2y) / denominator
    b = (p1x*p2y - p2x*p1y) / denominator
    y = a*x + b
    return x, y
  return (p1x, p2x), (p1y, p2y)  # matplotlib format


def cubic(point1, angle1, point2, angle2, samples=10):
  """Cubic line defined by two end points and the rotation in radians
  at the points.

  :param point1: one end point
  :type point1: tuple
  :param angle1: the slope at the one end point
  :type angle1: int or float
  :param point2: other end point
  :type point2: tuple
  :param angle2: the slope at the other end point
  :type angle2: int or float
  :param samples: number of sampling points (default = 10)
  :type samples: int

  :returns: ([sample_point1_x, sample_point2_x, ...],
    [sample_points1_y, sample_point2_y, ...])
  :rtype: tuple
  """
  p1x, p1y = point1
  p2x, p2y = point2

  x = np.linspace(p1x, p2x, num=samples)

  p1ys = math.tan(angle1)
  p2ys = math.tan(angle2)
  a = (p1x*p1ys + p1x*p2ys - p2x*p1ys - p2x*p2ys - 2*p1y + 2*p2y)/(p1x**3 - 3*p1x**2*p2x + 3*p1x*p2x**2 - p2x**3)
  b = (- p1x**2*p1ys - 2*p1x**2*p2ys - p1x*p2x*p1ys + p1x*p2x*p2ys + 3*p1x*p1y - 3*p1x*p2y + 2*p2x**2*p1ys + p2x**2*p2ys + 3*p2x*p1y - 3*p2x*p2y)/(p1x**3 - 3*p1x**2*p2x + 3*p1x*p2x**2 - p2x**3)
  c = (p1x**3*p2ys + 2*p1x**2*p2x*p1ys + p1x**2*p2x*p2ys - p1x*p2x**2*p1ys - 2*p1x*p2x**2*p2ys - 6*p1x*p2x*p1y + 6*p1x*p2x*p2y - p2x**3*p1ys)/(p1x**3 - 3*p1x**2*p2x + 3*p1x*p2x**2 - p2x**3)
  d = (- p1x**3*p2x*p2ys + p1x**3*p2y - p1x**2*p2x**2*p1ys + p1x**2*p2x**2*p2ys - 3*p1x**2*p2x*p2y + p1x*p2x**3*p1ys + 3*p1x*p2x**2*p1y - p2x**3*p1y)/(p1x**3 - 3*p1x**2*p2x + 3*p1x*p2x**2 - p2x**3)
  y = a*x**3 + b*x**2 + c*x + d
  return x, y


def cubic_deg(point1, angle1, point2, angle2):
  """Cubic line defined by two end points and the roation in degree
  at the points.

  :param point1: one end point
  :type point1: tuple
  :param angle1: the slope at the one end point
  :type angle1: int or float
  :param point2: other end point
  :type point2: tuple
  :param angle2: the slope at the other end point
  :type angle2: int or float

  :returns: ([sample_point1_x, sample_point2_x, ...],
    [sample_points1_y, sample_point2_y, ...])
  :rtype: tuple

  .. seealso::
    :meth:`cubic`
  """
  return cubic(point1, angle1 * math.pi/180, point2, angle2 * math.pi/180)