# -*- coding: utf-8 -*-
#
# Copyright © 2020 MSC.Software Corporation
#  
# Permission to use Scripts
#  
# Permission is granted to MSC.Software customers with valid licenses of the MSC Apex™ software 
# to use script files posted on this page (collectively “Scripts”) solely for customers’ own authorized 
# internal use in conjunction with licensed use of the MSC Apex software.  
# No other use is permitted.  No license to MSC Apex or other software is granted under this Permission.
#  
# SCRIPTS ARE PROVIDED ON AN "AS-IS" BASIS, WITH NO WARRANTIES OF ANY KIND.  USE OF THE SCRIPTS IS AT CUSTOMER'S SOLE RISK.
#  
# IN NO EVENT SHALL MSC.SOFTWARE CORPORATION, ITS AFFILIATES, SUBSIDIARIES, SUPPLIERS, AGENTS, REPRESENTATIVES OR EMPLOYEES BE LIABLE FOR ANY 
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL OR CONSEQUENTIAL DAMAGES HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, TORT, 
# STRICT LIABILITY, OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION USE OF THE SCRIPTS, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
#  
# Customer shall indemnify, defend and hold MSC.Software Corporation, its affiliates, subsidiaries, suppliers, agents, 
# employees or representatives (collectively, “Indemnified Parties”) harmless from and against any third-party claims 
# brought against Indemnified Parties rising out of or in connection with Customer’s use of the Scripts.
#
"""
Created on Thu Apr 21 12:55:09 2016
Miscellaneous routines for composite design explorer
@author: e1851
"""
try:
    import numpy
except:
    print("numpy not installed")
import apex
import math


def createABD(plyList, mat1Dict,mat8Dict):
# creates ABD matrix for given list of plies
    debug = False
    for ply in plyList:
        angle = ply.getAngle()
        material = (ply.getMaterial()).getName()
        thickness = ply.getThickness()
        # print(f"ply: {angle},{material},{thickness}")
    a = numpy.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])
    b = numpy.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])
    d = numpy.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])
# get total thickness
    t_total = 0.0
    for ply in plyList:
        t_total += ply.getThickness()
# set starting z at -t/2
    zbot= -t_total/2
    for ply in plyList:
        angle = ply.getAngle()
        thickness = ply.getThickness()
        # find CDE material object in dictionaries
        matId = (ply.getMaterial()).getName()
        # print(f"58ply: {angle},{material},{thickness}")
        material = mat1Dict.get(matId)
        if material is None:
            material = mat8Dict.get(matId)
        ztop = zbot + thickness
        q = material.q
        q_bar = transform_q(q, angle)
        a += q_bar * thickness
        b += 0.5 * q_bar * (ztop * ztop - zbot * zbot)
        d += q_bar / 3.0 * (ztop * ztop * ztop - zbot * zbot * zbot)
        zbot = ztop
    #  print(f"ABD {a} {b} {d}")
    return(a, b, d)


def create_polar_stiffness(k, plot_label="polar plot"):
    angle_list = []
    k_list = []
    theta = 0.0
    flex = numpy.linalg.inv(k)
    while theta <= 360.0:
        theta_rad = theta/90.0*math.acos(0.0)
        # apply unit force at angle theta
        f_vec_transpose = mohrs_rotate(1.0, 0.0, 0.0, theta_rad)
        f_vec = numpy.matrix(f_vec_transpose).transpose()
        # strain vector calculated
        strain = numpy.dot(flex, f_vec)
#        strain[2] = abs(strain[2])
        rotated_strain = mohrs_rotate(strain[0], strain[1], -strain[2]/2,
                                      -theta_rad)
        k_theta = float(1.0/rotated_strain[0])
#        print (theta, rotated_strain[0])
        angle_list.append(theta_rad)
        k_list.append(k_theta)
        if theta == 0.0:
            k_max = k_theta
            k_max_angle = 0.0
            k_min = k_theta
            k_min_angle = 0.0
        else:
            if k_theta > k_max:
                k_max = k_theta
                k_max_angle = theta
            if k_theta < k_min:
                k_min = k_theta
                k_min_angle = theta
        theta += 1
    stiff = (angle_list, k_list)
    print(f"{plot_label} Max stiffness occurs at {k_max_angle} degrees")
    print(f"{plot_label} Min stiffness occurs at {k_min_angle} degrees")
    print(f"{plot_label} max/min stiffness ratio = {k_max/k_min}")
    return stiff


def transform_q(q, theta):
    c = math.cos(theta/90.0*math.acos(0.0))
    c2 = c * c
    s = math.sin(theta/90.0*math.acos(0.0))
    s2 = s * s
    sc = s * c
    t = numpy.matrix([[c2, s2, sc], [s2, c2, -sc], [-2.0*sc, 2.0*sc, c2-s2]])
    t_trans = t.transpose()
    temp = numpy.dot(t_trans, q)
    qbar = numpy.dot(temp, t)
    return qbar


def mohrs_rotate(f1, f2, f12, angle):
    ''' uses mohr's circle to rotate forces, stresses or strains by angle.
    For strain ensure tensor shear strain is used.
  Angle is in radians'''
    f_ave = (f1 + f2) / 2
    radius = math.sqrt(f12 * f12 + (f1 - f_ave) ** 2)
    if radius <= 1.0e-6*f_ave:
        phi = 0.0
    else:
        phi = -math.atan2(f12, (f1 - f_ave))
    #print(f"phi: {phi}")
    angle2 = 2.0 * angle + phi
    fx = f_ave + radius * math.cos(angle2)
    fy = f_ave - radius * math.cos(angle2)
    fxy = radius * math.sin(angle2)
    return [fx, fy, fxy]
