'''
Purpose: Computing laminate properties, and transformation matrices
Author: Compiled from code by Nils Petter Vedvik
URL: https://folk.ntnu.no/nilspv/TMM4175
Requires: NumPy
'''

import numpy as np

def S2D(m):
    return np.array([[        1/m['E1'], -m['v12']/m['E1'],          0],
                     [-m['v12']/m['E1'],         1/m['E2'],          0],
                     [                0,                 0, 1/m['G12']]], float)
def Q2D(m):
    S=S2D(m)
    return np.linalg.inv(S)

def T2Ds(a):
    c,s = np.cos(np.radians(a)), np.sin(np.radians(a))
    return np.array([[ c*c ,  s*s ,   2*c*s],
                     [ s*s ,  c*c ,  -2*c*s],
                     [-c*s,   c*s , c*c-s*s]], float)

def T2De(a):
    c,s = np.cos(np.radians(a)), np.sin(np.radians(a))
    return np.array([[   c*c,   s*s,     c*s ],
                     [   s*s,   c*c,    -c*s ],
                     [-2*c*s, 2*c*s, c*c-s*s ]], float)

def laminateThickness(layup):
    return sum([layer['thi'] for layer in layup])

def Q2Dtransform(Q,a):
    return np.dot(np.linalg.inv( T2Ds(a) ) , np.dot(Q, T2De(a)) )

def computeD(layup):
    D = np.zeros((3, 3))
    h_bot = -laminateThickness(layup)/2
    for layer in layup:
        mat = layer['mat']
        t = layer['ori']
        h_top = h_bot + layer['thi']

        Q_prime = np.dot(np.linalg.inv(T2Ds(t)), np.dot(Q2D(mat), T2De(t)))
        D += (1 / 3) * Q_prime * (h_top**3 - h_bot**3)

        h_bot = h_top #bot of next layer is top of this layer

    return D


def computeA(layup):
    A=np.zeros((3,3),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q, a)
        htop = hbot + layer['thi']   # top of current layer
        A += Qt*(htop-hbot)
        hbot = htop                    # for the next layer
    return A