# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

from modsim import *

UNITS.Quantity.__init__


# +
class _Vector(UNITS.Quantity):
    """Represented as a Pint Quantity with a NumPy array
    
    x, y, z, mag, mag2, and angle are accessible as attributes.
    
    Supports vector operations hat, dot, cross, proj, and comp.
    """
    
    @property
    def x(self):
        """Returns the x component with units."""
        return self[0]

    @property
    def y(self):
        """Returns the y component with units."""
        return self[1]

    @property
    def z(self):
        """Returns the z component with units."""
        return self[2]

    @property
    def mag(self):
        """Returns the magnitude with units."""
        return np.sqrt(np.dot(self, self)) * self.units

    @property
    def mag2(self):
        """Returns the magnitude squared with units."""
        return np.dot(self, self) * self.units

    @property
    def angle(self):
        """Returns the angle between self and the positive x axis."""
        return np.arctan2(self.y, self.x)

    def polar(self):
        """Returns magnitude and angle."""
        return self.mag, self.angle

    def hat(self):
        """Returns the unit vector in the direction of self."""
        """
        """
        return self / self.mag * self.units

    def dot(self, other):
        """Returns the dot product of self and other."""
        """
        """
        return np.dot(self, other) * self.units * other.units

    def cross(self, other):
        """Returns the cross product of self and other."""
        """
        """
        return np.cross(self, other) * self.units * other.units

    def proj(self, other):
        """Returns the projection of self onto other."""
        """
        """
        return np.dot(self, other) * other.hat()

    def comp(self, other):
        """Returns the magnitude of the projection of self onto other."""
        """
        """
        return np.dot(self, other.hat()) * other.units

    def dist(self, other):
        """Euclidean distance from self to other, with units."""
        diff = self - other
        return diff.mag

    def diff_angle(self, other):
        """
        """
        #TODO: see http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/
        raise NotImplementedError()
        
        
def Vector(*args, units=None):
    # if there's only one argument, it should be iterable
    if len(args) == 1:
        args = args[0]
        
        # if it's a series, pull out the values
        if isinstance(args, Series):
            args = args.values
        
    # see if any of the arguments have unit; if so, save the first one
    for elt in args:
        found_units = getattr(elt, 'units', None)
        if found_units:
            break
            
    if found_units:
        # if there are units, remove them
        args = [getattr(elt, 'magnitude', elt) for elt in args]
    
    # if the units keyword is provided, it overrides the units in args
    if units is not None:
        found_units = units
    
    return _Vector(args, found_units)


# -

m = UNITS.m
N = UNITS.newton

a = Vector(3, 4)
a

a = Vector(3, 4, units=m)
a

a = Vector([3, 4])
a

a = Vector(3, 4*m)
a

a = Vector([3, 4*m])
a

a = Vector(3, 4*m, units=s) * m
a

a = Vector(Series([3, 4]))
a

a = Vector(3, 4)
a

a = Vector([3, 4]) * m
type(a)

a

a * m

a / m

a[0]

a.x

a.y

a.mag

a.mag2

type(a)

type(a.mag)

c = a * a.mag
c

c = a / a.mag
c

c.x

ah = a.hat()
ah



print(ah)

ah.units





b = Vector([1, 2])

a.dot(b)

# +
N = UNITS.newton

b = Vector([1, 2]) * N
# -

a.dot(b)

c = a.cross(b)
c

c = a.proj(b)

c

a.comp(b)

x_hat = Vector([1, 0])
y_hat = Vector([0, 1])

a.comp(x_hat)

a.comp(y_hat)

a.angle

b.angle

Vector([1, 1]).angle / pi

Vector([-1, 1]).angle / pi

Vector([-1, -1]).angle / pi

Vector([1, -1]).angle / pi

Vector([1, -1]).polar()

x = [5, 3.5355, 0, -10]
x

y = [0, 3.5355, 10, 0]
y

theta, rho = cart2pol(x,y)

theta

rho

theta = [0, pi/4, pi/2, pi]

rho = [5, 5, 10, 10]

x2, y2 = pol2cart(theta,rho)

print(x2)

print(y2)

degree = UNITS.degree
kg = UNITS.kg
m = UNITS.m
s = UNITS.s

condition = Condition(x = 0 * m, 
                      y = 0 * m,
                      g = 9.8 * m/s**2,
                      mass = 145e-3 * kg,
                      diameter = 73e-3 * m,
                      rho = 1.2 * kg/m**3,
                      C_d = 0.3,
                      angle = 45 * degree,
                      velocity = 40 * m / s,
                      duration = 5 * s)


def make_system(condition):
    unpack(condition)
    
    theta = np.deg2rad(angle)
    vx, vy = pol2cart(theta, velocity)
    init = State(x=x, y=y, vx=vx, vy=vy)
    v = State(vx=vx, vy=vy)
    
    area = np.pi * (diameter/2)**2
    ts = linspace(0, duration, 101)
    
    return System(init=init, g=g, mass=mass, v=v,
                  area=area, rho=rho, C_d=C_d, ts=ts)


system = make_system(condition)
system.init

v = system.v.values
v





def slope_func(state, t, system):
    x, y, vx, vy = state
    unpack(system)
    
    a_grav = Vector(0, -g)

    v = Vector(vx, vy)
    
    f_drag = -rho * v.mag * v * C_d * area / 2
    a_drag = f_drag / mass

    a = a_grav + a_drag
    
    return v, a


slope_func(system.init, 0, system)

a = 5 * m

b = Vector([1, 2]) * s

b * a

a * b













StarcoderdataJupyterScriptsDedupFiltered