Pane.py

# -*- coding: utf-8 -*-
"""
This file is a part of mei_comp_2020

This file creates the triangular panes that make the object. 
Three (x, y, z) locations of the three vertices of the panes are 
needed in a 2D list. The order of points needs to be such that the 
cross product of  (point[1] – point[0]) and (point[2] – point[0]) 
is away from the center of the object. A density, thickness, and dt 
is also needed.

Requires the use of numpy.
"""
#==============================================================================
import numpy as np

class Pane(object):
    def __init__(self, location, density, thickness, dt, verbose = 0):
        self.points = np.array(location)
        self.den = density
        self.thick = thickness
        self.dt = dt
        self.verbose = verbose
        
        self.area = self.calcArea()
        self.mass = self.calcMass()
        self.CM = self.calcCM()
        self.nHat = self.calcNHat()
        
#==============================================================================

    """These definitions are used within this file's init values 
    to calculate needed values"""

    def calcArea(self):
        return(.5 * self.mag(np.cross((self.points[1] - self.points[0]), (self.points[2] - self.points[0]))))
        
    def calcMass(self):
        return(self.area * self.thick *self.den)

    def calcCM(self):
        return np.array(((self.points[0,0] + self.points[1,0] + self.points[2,0])/3,
                         (self.points[0,1] + self.points[1,1] + self.points[2,1])/3,
                         (self.points[0,2] + self.points[1,2] + self.points[2,2])/3))
 
    def calcNHat(self):
        n = np.cross((self.points[1] - self.points[0]), (self.points[2] - self.points[0]))
        magN = np.sqrt(np.square(n[0]) + np.square(n[1]) + np.square(n[2]))
        return(n / magN)
        
#==============================================================================
        
    """this definition is called by Object.py to calculate the force on the pane"""

    def calcForceResistive(self, vel, denAir, mAir, theta):
        nHat = self.rotate(theta, self.nHat)

        velMag = np.sqrt(np.square(vel[0]) + np.square(vel[1]) + np.square(vel[2]))
        #think about this small number
        if (velMag <= 1e-8):
            velHat = np.array((0., 0., 0.))
        else:    
            velHat = vel / velMag
        dot = np.dot(velHat, nHat)
        
        if (self.verbose >= 3):
            print("In calcForceResistive, vel=[{:.3g}, {:.3g}, {:.3g}], "
                  "nHat=[{:6.3f}, {:6.3f}, {:6.3f}]"
                  .format(vel[0], vel[1], vel[2], nHat[0], nHat[1], nHat[2]))
        
        if dot <= 0:
            return np.zeros(3)
        
        useableArea = (self.area * dot)
        numParticles = useableArea * (self.dt * velMag) * denAir
        dpp = -2 * velMag * np.dot(velHat, nHat) * nHat * mAir * numParticles

        return (dpp / self.dt)
    
#==============================================================================

    """These are definitions to be used within this file for ease of use"""
    
    def rotate(self, theta, xy):
        return (np.array((((xy[0] * np.cos(theta)) - (xy[1] * np.sin(theta))), 
                         ((xy[0] * np.sin(theta)) + (xy[1] * np.cos(theta))), 0)))
    
    def mag(self, r):
        R = 0
        for i in range(len(r)):
            R += np.square(r[i])
        return np.sqrt(R)