From cd249eae043bd51cd7c91359f0a22dad22bde14d Mon Sep 17 00:00:00 2001 From: Aleksandr Slatov <115647936+alesanter@users.noreply.github.com> Date: Fri, 12 Jul 2024 02:50:20 +0400 Subject: [PATCH] Law: displacement of falling body (#845) * Add convenience function for vectors * Create law outline and test * Refactor vector addition and subtraction * Add law description * Fix vector arithmetic functions * Fix vector addition function usage * Add note about Big O --- symplyphysics/core/vectors/arithmetics.py | 71 +++++++----- .../vector/falling_body_displacement.py | 103 ++++++++++++++++++ test/core/vectors/arithmetics_test.py | 2 +- .../vector/falling_body_displacement_test.py | 88 +++++++++++++++ 4 files changed, 236 insertions(+), 28 deletions(-) create mode 100644 symplyphysics/laws/gravity/vector/falling_body_displacement.py create mode 100644 test/gravity/vector/falling_body_displacement_test.py diff --git a/symplyphysics/core/vectors/arithmetics.py b/symplyphysics/core/vectors/arithmetics.py index 7fe40e6c1..3fc8fd885 100644 --- a/symplyphysics/core/vectors/arithmetics.py +++ b/symplyphysics/core/vectors/arithmetics.py @@ -38,36 +38,53 @@ def equal_vectors(vector_left: Vector, vector_right: Vector) -> bool: return True -# Sum of two vectors -# Sum of two vectors can be seen as a diagonal of the parallelogram, where vectors are adjacent sides of this parallelogram. -# To subtract vectors, multiply one of the vectors to -1 and add them. +# Add any number of vectors #NOTE: adding two non-cartesian vectors is not a trivial task. We suggest to convert them to cartesian # vectors, add them and convert back. -def add_cartesian_vectors(vector_left: Vector, vector_right: Vector) -> Vector: - if vector_left.coordinate_system != vector_right.coordinate_system: - raise TypeError( - f"Different coordinate systems in vectors: {str(vector_left.coordinate_system)} vs {str(vector_right.coordinate_system)}" - ) - if vector_left.coordinate_system.coord_system_type != CoordinateSystem.System.CARTESIAN: - coord_name_from = CoordinateSystem.system_to_transformation_name( - vector_left.coordinate_system.coord_system_type) - raise ValueError( - f"Addition is only supported for cartesian coordinates: got {coord_name_from}") - if vector_right.coordinate_system.coord_system_type != CoordinateSystem.System.CARTESIAN: - coord_name_from = CoordinateSystem.system_to_transformation_name( - vector_right.coordinate_system.coord_system_type) - raise ValueError( - f"Addition is only supported for cartesian coordinates: got {coord_name_from}") - (list_left_extended, list_right_extended) = _extend_two_vectors(vector_left, vector_right) - result = list( - map(lambda lr: sympify(lr[0] + lr[1]), zip(list_left_extended, list_right_extended))) - return Vector(result, vector_left.coordinate_system) - - -def subtract_cartesian_vectors(vector_left: Vector, vector_right: Vector) -> Vector: +def add_cartesian_vectors(*vectors: Vector) -> Vector: + def add_two_cartesian_vectors(vector_left: Vector, vector_right: Vector) -> Vector: + if vector_left.coordinate_system != vector_right.coordinate_system: + raise ValueError( + "Vectors must have the same coordinate system, " + f"got {vector_left.coordinate_system} and {vector_right.coordinate_system}" + ) + + left_type = vector_left.coordinate_system.coord_system_type + right_type = vector_right.coordinate_system.coord_system_type + + if left_type != CoordinateSystem.System.CARTESIAN or right_type != CoordinateSystem.System.CARTESIAN: + left_name = CoordinateSystem.system_to_transformation_name(left_type) + right_name = CoordinateSystem.system_to_transformation_name(right_type) + raise ValueError(f"Expect two Cartesian vectors, got {left_name} and {right_name}") + + (list_left_extended, list_right_extended) = _extend_two_vectors(vector_left, vector_right) + result = [ + sympify(left_component + right_component) + for (left_component, right_component) + in zip(list_left_extended, list_right_extended) + ] + + return Vector(result, vector_left.coordinate_system) + + if len(vectors) < 1: + raise ValueError("Provide at least one vector") + + result = vectors[0] + for vector in vectors[1:]: + result = add_two_cartesian_vectors(result, vector) + return result + + +# Subtract vectors from the first given vector +def subtract_cartesian_vectors(*vectors: Vector) -> Vector: + if len(vectors) < 2: + raise ValueError("Provide at least two vectors") + + vector_minuend = vectors[0] + vector_subtrahend = add_cartesian_vectors(*vectors[1:]) return add_cartesian_vectors( - vector_left, - scale_vector(-1, vector_right) + vector_minuend, + scale_vector(-1, vector_subtrahend), ) diff --git a/symplyphysics/laws/gravity/vector/falling_body_displacement.py b/symplyphysics/laws/gravity/vector/falling_body_displacement.py new file mode 100644 index 000000000..702d909ab --- /dev/null +++ b/symplyphysics/laws/gravity/vector/falling_body_displacement.py @@ -0,0 +1,103 @@ +from sympy import Expr, Rational +from symplyphysics import ( + angle_type, + units, + Quantity, + Vector, + QuantityVector, + validate_input, + validate_output, + scale_vector, + add_cartesian_vectors, + cross_cartesian_vectors, +) + +# Description +## Suppose a reference frame S' is fixed to a moving object A (e.g., Earth) and some body B moving freely +## (i.e. the sum of external non-gravitational forces acting on it is zero). In the case of an inertial +## frame of reference, the displacement of body B from the starting position would follow the usual rule +## `s = v0 * t + (g / 2) * t**2`. But in the case of non-inertial frames, we have to take the Coriolis and +## the centrifugal force into account as well, which results into the following series: + +# Conditions +## - The sum `F` of all other, non-gravitational forces acting on body B is 0. +## - `g` is independent of coordinates. + +# Notes +## - Note that the series is truncated at the fifth term. More terms can be obtained by plugging the result +## into the equation of motion `a = g + cross(v, w)` and integrating it over time. + +# Law: s = v0 * t +# + t**2 * (g / 2 + cross(v0, w)) +# + t**3 / 3 * (cross(g, w) + 2 * cross(cross(v0, w), w)) +# + t**4 / 6 * cross(cross(g, w), w) +# + O(t**5) +## t - time +## s - vector of displacement of body B in S' +## v0 - vector of initial velocity of body B +## w - pseudovector of angular velocity of rotation of moving frame S' +## g - vector of acceleration due to gravity (or free fall acceleration) of body B +## cross(a, b) - cross product between vectors `a` and `b` +## O - [big O](https://en.wikipedia.org/wiki/Big_O_notation) as per the formal mathematical definition +## in the limit of t approaching infinity + + +def displacement_law( + time_: Expr, + initial_velocity_: Vector, + angular_velocity_: Vector, + acceleration_due_to_gravity_: Vector, +) -> Vector: + initial_cross_angular = cross_cartesian_vectors( + initial_velocity_, + angular_velocity_, + ) + + acceleration_cross_angular = cross_cartesian_vectors( + acceleration_due_to_gravity_, + angular_velocity_, + ) + + return add_cartesian_vectors( + scale_vector(time_, initial_velocity_), + scale_vector( + time_**2, + add_cartesian_vectors( + scale_vector(Rational(1, 2), acceleration_due_to_gravity_), + initial_cross_angular, + ), + ), + scale_vector( + time_**3 / 3, + add_cartesian_vectors( + acceleration_cross_angular, + scale_vector(2, cross_cartesian_vectors(initial_cross_angular, angular_velocity_)), + ), + ), + scale_vector( + time_**4 / 6, + cross_cartesian_vectors(acceleration_cross_angular, angular_velocity_) + ) + ) + + +@validate_input( + time_=units.time, + initial_velocity_=units.velocity, + angular_velocity_=angle_type / units.time, + acceleration_due_to_gravity_=units.acceleration, +) +@validate_output(units.length) +def calculate_displacement( + time_: Quantity, + initial_velocity_: QuantityVector, + angular_velocity_: QuantityVector, + acceleration_due_to_gravity_: QuantityVector, +) -> QuantityVector: + result_vector = displacement_law( + time_, + initial_velocity_.to_base_vector(), + angular_velocity_.to_base_vector(), + acceleration_due_to_gravity_.to_base_vector(), + ) + return QuantityVector.from_base_vector(result_vector) diff --git a/test/core/vectors/arithmetics_test.py b/test/core/vectors/arithmetics_test.py index 6020623e4..a58bf01ce 100644 --- a/test/core/vectors/arithmetics_test.py +++ b/test/core/vectors/arithmetics_test.py @@ -92,7 +92,7 @@ def test_rebased_add_vectors(test_args: Args) -> None: def test_invalid_add_vectors(test_args: Args) -> None: C1 = CoordinateSystem() - with raises(TypeError): + with raises(ValueError): add_cartesian_vectors(Vector([1, 2], test_args.C), Vector([1, 2], C1)) # non-cartesian addition is not supported C2 = CoordinateSystem(CoordinateSystem.System.CYLINDRICAL) diff --git a/test/gravity/vector/falling_body_displacement_test.py b/test/gravity/vector/falling_body_displacement_test.py new file mode 100644 index 000000000..02c7dee50 --- /dev/null +++ b/test/gravity/vector/falling_body_displacement_test.py @@ -0,0 +1,88 @@ +from collections import namedtuple +from pytest import fixture, raises +from symplyphysics import ( + errors, + units, + Quantity, + QuantityVector, + assert_equal_vectors, +) +from symplyphysics.laws.gravity.vector import falling_body_displacement as law + +Args = namedtuple("Args", "t v0 w g") + + +@fixture(name="test_args") +def test_args_fixture() -> Args: + t = Quantity(5 * units.second) + + v_unit = units.meter / units.second + v0 = QuantityVector([v_unit, v_unit, v_unit]) + + w = QuantityVector([0, -7.27e-5 * units.radian / units.second, 0]) + + g = QuantityVector([0, 0, -9.81 * units.meter / units.second**2]) + + return Args(t=t, v0=v0, w=w, g=g) + + +def test_law(test_args: Args) -> None: + result = law.calculate_displacement(test_args.t, test_args.v0, test_args.w, test_args.g) + assert_equal_vectors( + result, + QuantityVector([4.97 * units.meter, 5.00 * units.meter, -118.0 * units.meter]), + tolerance=4e-3, + ) + + +def test_bad_time(test_args: Args) -> None: + tb = Quantity(units.coulomb) + with raises(errors.UnitsError): + law.calculate_displacement(tb, test_args.v0, test_args.w, test_args.g) + with raises(TypeError): + law.calculate_displacement(100, test_args.v0, test_args.w, test_args.g) + + +def test_bad_linear_velocity(test_args: Args) -> None: + vb_vector = QuantityVector([units.coulomb]) + with raises(errors.UnitsError): + law.calculate_displacement(test_args.t, vb_vector, test_args.w, test_args.g) + + vb_scalar = Quantity(units.meter / units.second) + with raises(AttributeError): + law.calculate_displacement(test_args.t, vb_scalar, test_args.w, test_args.g) + + with raises(TypeError): + law.calculate_displacement(test_args.t, 100, test_args.w, test_args.g) + with raises(TypeError): + law.calculate_displacement(test_args.t, [100], test_args.w, test_args.g) + + +def test_bad_angular_velocity(test_args: Args) -> None: + wb_vector = QuantityVector([units.coulomb]) + with raises(errors.UnitsError): + law.calculate_displacement(test_args.t, test_args.v0, wb_vector, test_args.g) + + wb_scalar = Quantity(units.radian / units.second) + with raises(AttributeError): + law.calculate_displacement(test_args.t, test_args.v0, wb_scalar, test_args.g) + + with raises(TypeError): + law.calculate_displacement(test_args.t, test_args.v0, 100, test_args.g) + with raises(TypeError): + law.calculate_displacement(test_args.t, test_args.v0, [100], test_args.g) + + +def test_bad_acceleration(test_args: Args) -> None: + gb_vector = QuantityVector([units.coulomb]) + with raises(errors.UnitsError): + law.calculate_displacement(test_args.t, test_args.v0, test_args.w, gb_vector) + + gb_scalar = Quantity(units.meter / units.second**2) + with raises(AttributeError): + law.calculate_displacement(test_args.t, test_args.v0, test_args.w, gb_scalar) + + with raises(TypeError): + law.calculate_displacement(test_args.t, test_args.v0, test_args.w, 100) + with raises(TypeError): + law.calculate_displacement(test_args.t, test_args.v0, test_args.w, [100])