cohete3/physics.py at master · kaz2700/cohete3 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
from  coords_math import radial_to_cartesian, cartesian_to_radial, norm, get_normal_vector
from physics_constants import G, EARTH_MASS, ATMOSPHERIC_PRESSURE_EARTH_SURFACE, EARTH_ANGULAR_VELOCITY, WIND_VELOCITY
import numpy as np
from parameters import dimensions

#physics related operations

#sums all forces and divides by the mass
def get_acceleration(rocket, add_thrust):
    forces = np.array([0, 0, 0])
    forces = forces + np.array(get_gravity_force(rocket))
    forces = forces + np.array(get_friction_force(rocket))
    forces = forces + np.array(get_coriolis_force(rocket))
    forces = forces + np.array(get_centrifuge_force(rocket))
    if add_thrust:
        forces = forces + np.array(get_thrust_force(rocket))
    return forces / rocket.get_mass()

#gets the rocket from the thrust and makes the force push in the direction of the vector pointing to the rocket (perpendicular to the surface of the earth)
def get_thrust_force(rocket):
    thrust_force = [0, 0, 0]
    rocket_thrust = rocket.get_thrust()
    normal_vector = get_normal_vector(rocket.positions[-1])

    for dim in range(dimensions):
        thrust_force += normal_vector[dim] * rocket_thrust

    return thrust_force

#gets the gravity force
def get_gravity_force(rocket):
    position = rocket.positions[-1]
    position_spherical_coords = cartesian_to_radial(position[0], position[1], position[2])
    gravity_force = -G * EARTH_MASS * rocket.get_mass() / rocket.get_radial_distance()**2

    return radial_to_cartesian(gravity_force, position_spherical_coords[1], position_spherical_coords[2]) #todo landa and alpha of rocket

#get atmospheric pressure by the formula provided
def get_atmospheric_pressure(rocket):
    if rocket.get_height() > 44000:
        return 0

    g = 9.80665 #m s-2
    kappa = 0.0034836177811575927 #m2 K-1 s-2
    L_B = 0.0065 #K m-1
    T_0 = 288.15 #K
    exponent = 1 - g*kappa/L_B

    pressure = ATMOSPHERIC_PRESSURE_EARTH_SURFACE * (T_0 / (T_0 - L_B*rocket.get_height()))**exponent

    return pressure

#get the friction force using the atmospheric pressure and with velocity relative to the wind
def get_friction_force(rocket):
    friction_force = [0, 0, 0]

    position_rocket = rocket.positions[-1]
    velocity_rocket = rocket.velocities[-1]

    velocity_wind = WIND_VELOCITY(position_rocket)

    velocity = [0, 0, 0]
    for dimension in range(dimensions):
        velocity[dimension] += velocity_rocket[dimension] - velocity_wind[dimension]

    velocity_norm = norm(velocity)

    area = 0
    for stage in rocket.stages:
        if not stage.is_detached: #if its detached it wont cause friction
            area += stage.Drag_Area * stage.num_of_boosters

    for dimension in range(dimensions):
        friction_force[dimension] += 0.5 * get_atmospheric_pressure(rocket) * area * -velocity[dimension] * velocity_norm

    return friction_force

#both the coriolis and the centrifuge forces are inertial forces caused because of the rotation of the system
def get_coriolis_force(rocket):
    coriolis_force = [0, 0, 0]

    velocity = rocket.velocities[-1]
    mass = rocket.get_mass()

    #with cross product
    coriolis_force[0] = -2 * mass * EARTH_ANGULAR_VELOCITY * velocity[1]
    coriolis_force[1] = -2 * mass * EARTH_ANGULAR_VELOCITY * velocity[0]

    return coriolis_force

def get_centrifuge_force(rocket):
    centrifuge_force = [0, 0, 0]
    position = rocket.positions[-1]
    mass = rocket.get_mass()

    centrifuge_force[0] = -mass * EARTH_ANGULAR_VELOCITY ** 2 * position[0]
    centrifuge_force[1] = -mass * EARTH_ANGULAR_VELOCITY ** 2 * position[1]

    return centrifuge_force