#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np


def remove_nan(xyzs): return xyzs[~np.isnan(xyzs).any(axis = 1)]


def measure_twocores(core_xyz_ref, core_xyz_tar):
    ''' Measure the following aspects of two helical cores.
        - Interhelical distance vector between the centers.
        - Interhelical angle (0-90 degree)
    '''
    # Obtain the centers...
    center_ref = np.nanmean(core_xyz_ref, axis = 0)
    center_tar = np.nanmean(core_xyz_tar, axis = 0)

    # Construct the interhelical distance vector...
    ih_dvec = center_tar - center_ref

    # Calculate the length of interhelical distance vector...
    norm_ih_dvec = np.linalg.norm(ih_dvec)

    # Obtain the helical core vectors...
    core_xyz_ref_nonan = remove_nan(core_xyz_ref)
    core_xyz_tar_nonan = remove_nan(core_xyz_tar)
    core_vec_ref = core_xyz_ref_nonan[-1] - core_xyz_ref_nonan[0]
    core_vec_tar = core_xyz_tar_nonan[-1] - core_xyz_tar_nonan[0]

    # Calculate the interhelical angle...
    core_vec_ref_unit = core_vec_ref / np.linalg.norm(core_vec_ref)
    core_vec_tar_unit = core_vec_tar / np.linalg.norm(core_vec_tar)
    ih_ang = np.arccos( np.dot(core_vec_ref_unit, core_vec_tar_unit) )

    return ih_dvec, norm_ih_dvec, core_vec_ref_unit, core_vec_tar_unit, ih_ang


def calc_interangle(core_xyz_ref, core_xyz_tar):
    ''' Measure the following aspects of two helical cores.
        - Interhelical angle (0-90 degree)
    '''
    # Obtain the helical core vectors...
    core_xyz_ref_nonan = remove_nan(core_xyz_ref)
    core_xyz_tar_nonan = remove_nan(core_xyz_tar)
    core_vec_ref = core_xyz_ref_nonan[-1] - core_xyz_ref_nonan[0]
    core_vec_tar = core_xyz_tar_nonan[-1] - core_xyz_tar_nonan[0]

    # Calculate the interhelical angle...
    core_vec_ref_unit = core_vec_ref / np.linalg.norm(core_vec_ref)
    core_vec_tar_unit = core_vec_tar / np.linalg.norm(core_vec_tar)
    inter_angle = np.arccos( np.dot(core_vec_ref_unit, core_vec_tar_unit) )

    if inter_angle > np.pi / 2.0: inter_angle = np.pi - inter_angle

    return inter_angle


def calc_interdist(core_xyz_ref, core_xyz_tar):
    ''' Measure the following aspects of two helical cores.
        - Interhelical distance vector between the centers.

        Refers to http://geomalgorithms.com/a07-_distance.html for the method.
        Q is ref, P is tar.  
    '''
    # Obtain the helical core vectors...
    core_xyz_ref_nonan = remove_nan(core_xyz_ref)
    core_xyz_tar_nonan = remove_nan(core_xyz_tar)
    core_vec_ref = core_xyz_ref_nonan[-1] - core_xyz_ref_nonan[0]
    core_vec_tar = core_xyz_tar_nonan[-1] - core_xyz_tar_nonan[0]

    # Obtain the starting point...
    q0 = core_xyz_ref_nonan[0]
    p0 = core_xyz_tar_nonan[0]
    w0 = p0 - q0

    # Obtain the directional vector with magnitude...
    v = core_vec_ref
    u = core_vec_tar

    # Math part...
    a = np.dot(u, u)
    b = np.dot(u, v)
    c = np.dot(v, v)
    d = np.dot(u, w0)
    e = np.dot(v, w0)
    de = a * c - b * b    # Denominator
    if de == 0: sc, tc = 0, d / b
    else:       sc, tc = (b * e - c * d) / de, (a * e - b * d) / de

    # Calculate distance...
    wc = w0 + sc * u - tc * v
    inter_dist = np.linalg.norm(wc)

    return inter_dist

[((16, 17, 16, 51), 'numpy.nanmean', 'np.nanmean', (), '', True, 'import numpy as np\n'), ((17, 17, 17, 51), 'numpy.nanmean', 'np.nanmean', (), '', True, 'import numpy as np\n'), ((23, 19, 23, 42), 'numpy.linalg.norm', 'np.linalg.norm', ({(23, 34, 23, 41): 'ih_dvec'}, {}), '(ih_dvec)', True, 'import numpy as np\n'), ((82, 8, 82, 20), 'numpy.dot', 'np.dot', ({(82, 15, 82, 16): 'u', (82, 18, 82, 19): 'u'}, {}), '(u, u)', True, 'import numpy as np\n'), ((83, 8, 83, 20), 'numpy.dot', 'np.dot', ({(83, 15, 83, 16): 'u', (83, 18, 83, 19): 'v'}, {}), '(u, v)', True, 'import numpy as np\n'), ((84, 8, 84, 20), 'numpy.dot', 'np.dot', ({(84, 15, 84, 16): 'v', (84, 18, 84, 19): 'v'}, {}), '(v, v)', True, 'import numpy as np\n'), ((85, 8, 85, 21), 'numpy.dot', 'np.dot', ({(85, 15, 85, 16): 'u', (85, 18, 85, 20): 'w0'}, {}), '(u, w0)', True, 'import numpy as np\n'), ((86, 8, 86, 21), 'numpy.dot', 'np.dot', ({(86, 15, 86, 16): 'v', (86, 18, 86, 20): 'w0'}, {}), '(v, w0)', True, 'import numpy as np\n'), ((93, 17, 93, 35), 'numpy.linalg.norm', 'np.linalg.norm', ({(93, 32, 93, 34): 'wc'}, {}), '(wc)', True, 'import numpy as np\n'), ((32, 39, 32, 67), 'numpy.linalg.norm', 'np.linalg.norm', ({(32, 54, 32, 66): 'core_vec_ref'}, {}), '(core_vec_ref)', True, 'import numpy as np\n'), ((33, 39, 33, 67), 'numpy.linalg.norm', 'np.linalg.norm', ({(33, 54, 33, 66): 'core_vec_tar'}, {}), '(core_vec_tar)', True, 'import numpy as np\n'), ((34, 24, 34, 68), 'numpy.dot', 'np.dot', ({(34, 31, 34, 48): 'core_vec_ref_unit', (34, 50, 34, 67): 'core_vec_tar_unit'}, {}), '(core_vec_ref_unit, core_vec_tar_unit)', True, 'import numpy as np\n'), ((50, 39, 50, 67), 'numpy.linalg.norm', 'np.linalg.norm', ({(50, 54, 50, 66): 'core_vec_ref'}, {}), '(core_vec_ref)', True, 'import numpy as np\n'), ((51, 39, 51, 67), 'numpy.linalg.norm', 'np.linalg.norm', ({(51, 54, 51, 66): 'core_vec_tar'}, {}), '(core_vec_tar)', True, 'import numpy as np\n'), ((52, 29, 52, 73), 'numpy.dot', 'np.dot', ({(52, 36, 52, 53): 'core_vec_ref_unit', (52, 55, 52, 72): 'core_vec_tar_unit'}, {}), '(core_vec_ref_unit, core_vec_tar_unit)', True, 'import numpy as np\n'), ((7, 35, 7, 49), 'numpy.isnan', 'np.isnan', ({(7, 44, 7, 48): 'xyzs'}, {}), '(xyzs)', True, 'import numpy as np\n')]