hw10

2.py

import numpy as np


def dist(angle, m_prop):
    t = np.linspace(0, 10, 1001)
    vx = np.zeros([1001])
    vy = np.zeros([1001])
    x = np.zeros([1001])
    y = np.zeros([1001])

    A = 0.8
    g = -9.8
    initial_height = 5
    initial_velocity = 1500 * (m_prop / 65)**.45
    mass = 65
    C = 1.4
    rho = 1.225
    radians = angle * 2 * np.pi / 360

    y[0] = initial_height
    vx[0] = initial_velocity * np.cos(radians)
    vy[0] = initial_velocity * np.sin(radians)

    for j in range(1, len(t)):
        if y[j - 1] <= 0.0:
            vx[j] = vy[j] = vx[j - 1] = vy[
                j - 1] = y[j] = y[j - 1] = 0.0
            x[j - 1] = x[j - 2]
            x[j] = x[j - 1]
            continue

        v = np.sqrt(vx[j - 1]**2 + vy[j - 1]**2)
        ax = -(0.5 * rho * C * A / mass) * v**2 * (vx[j - 1] / v)
        ay = g - (0.5 * rho * C * A / mass) * v**2 * (vy[j - 1] / v)
        dt = t[1] - t[0]

        vy[j] = vy[j - 1] + ay * dt
        vx[j] = vx[j - 1] + ax * dt

        x[j] = x[j - 1] + vx[j] * dt
        y[j] = y[j - 1] + vy[j] * dt

    maxRow = x.max(0)
    return maxRow


print(dist(30, 0.10))

3.py

import numpy as np


def dist(angle, m_prop):
    t = np.linspace(0, 10, 1001)
    vx = np.zeros([1001])
    vy = np.zeros([1001])
    x = np.zeros([1001])
    y = np.zeros([1001])

    A = 0.8
    g = -9.8
    initial_height = 5
    initial_velocity = 1500 * (m_prop / 65)**.45
    mass = 65
    C = 1.4
    rho = 1.225
    radians = angle * 2 * np.pi / 360

    y[0] = initial_height
    vx[0] = initial_velocity * np.cos(radians)
    vy[0] = initial_velocity * np.sin(radians)

    for j in range(1, len(t)):
        if y[j - 1] <= 0.0:
            vx[j] = vy[j] = vx[j - 1] = vy[
                j - 1] = y[j] = y[j - 1] = 0.0
            x[j - 1] = x[j - 2]
            x[j] = x[j - 1]
            continue

        v = np.sqrt(vx[j - 1]**2 + vy[j - 1]**2)
        ax = -(0.5 * rho * C * A / mass) * v**2 * (vx[j - 1] / v)
        ay = g - (0.5 * rho * C * A / mass) * v**2 * (vy[j - 1] / v)
        dt = t[1] - t[0]

        vy[j] = vy[j - 1] + ay * dt
        vx[j] = vx[j - 1] + ax * dt

        x[j] = x[j - 1] + vx[j] * dt
        y[j] = y[j - 1] + vy[j] * dt

    maxRow = x.max(0)
    return maxRow


# print(dist(30, 0.10))

guess_dist   = 0    # m
guess_angle  = 0    # deg
guess_m_prop = 0.05 # kg

max_dist = 0.0
max_angle = 0
max_m_prop = 0

# change angle iteratively until improvement stops
old_guess_dist = guess_dist - 1
while guess_dist >= old_guess_dist:
    old_guess_dist = guess_dist
    guess_angle += 1
    guess_dist = dist( guess_angle,guess_m_prop )
    max_angle = guess_angle

# then change m_prop iteratively until improvement stops
old_guess_dist = guess_dist - 1
while guess_dist >= old_guess_dist:
    max_m_prop = guess_m_prop
    old_guess_dist = guess_dist
    guess_m_prop += 0.01
    if guess_m_prop > 0.12: break
    guess_dist = dist( max_angle,guess_m_prop )
    print(guess_m_prop, guess_dist)
max_dist = guess_dist

# now you should know max_dist, max_angle, max_m_prop

4.py

from scipy.optimize import minimize
import math
def f_r(x):
    return 20 + x**2 - 0.1 * math.cos(math.pi * 2 * x)
r = minimize(f_r, 1)
xstar = r['x'][0]
fstar = r['fun']