How to plot the motion of a projectile under the effect of gravity, buoyancy and air resistance?












1















I am trying to make a plot of a projectile motion of a mass which is under the effect of gravitational, buoyancy and drag force. Basically, I want to show effects of the buoyancy and drag force on flight distance, flight time and velocity change on plotting.



import matplotlib.pyplot as plt

import numpy as np



V_initial = 30 # m/s

theta = np.pi/6 # 30

g = 3.711

m =1

C = 0.47

r = 0.5

S = np.pi*pow(r, 2)

ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)

t = np.linspace(0, t_flight, 200)



# Drag force

Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)



# Buoyant Force

Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))



x_loc = 

y_loc = 



for time in t:

    x = V_initial*time*np.cos(theta)

        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)

    x_loc.append(x)

    y_loc.append(y)



x_vel = 

y_vel = 

for time in t:

    vx = V_initial*np.cos(theta)

    vy = V_initial*np.sin(theta) - g*time

    x_vel.append(vx)

    y_vel.append(vy)





v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]



tau = 

for velocity in v_ch:

    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)

        tau.append(Ft)



buoy = 

for velocity in v_ch:

    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

    buoy.append(Fb)


after this point, I couldn't figure out how to plot to a projectile motion under this forces. In other words, I'm trying to compare the projectile motion of the mass under three circumstances




  1. Mass only under the effect of gravity

  2. Mass under the effect of gravity and air resistance

  3. Mass under the effect of gravity, air resistance, and buoyancy






python numpy matplotlib physics







share|improve this question




















  • 1





    You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

    – Bazingaa
    Jan 2 at 8:19













  • From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

    – Mstaino
    Jan 2 at 13:05











  • I edited my question and tried to explain more clearly.

    – Birkan Emrem
    Jan 2 at 13:28
















1















I am trying to make a plot of a projectile motion of a mass which is under the effect of gravitational, buoyancy and drag force. Basically, I want to show effects of the buoyancy and drag force on flight distance, flight time and velocity change on plotting.



import matplotlib.pyplot as plt

import numpy as np



V_initial = 30 # m/s

theta = np.pi/6 # 30

g = 3.711

m =1

C = 0.47

r = 0.5

S = np.pi*pow(r, 2)

ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)

t = np.linspace(0, t_flight, 200)



# Drag force

Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)



# Buoyant Force

Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))



x_loc = 

y_loc = 



for time in t:

    x = V_initial*time*np.cos(theta)

        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)

    x_loc.append(x)

    y_loc.append(y)



x_vel = 

y_vel = 

for time in t:

    vx = V_initial*np.cos(theta)

    vy = V_initial*np.sin(theta) - g*time

    x_vel.append(vx)

    y_vel.append(vy)





v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]



tau = 

for velocity in v_ch:

    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)

        tau.append(Ft)



buoy = 

for velocity in v_ch:

    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

    buoy.append(Fb)


after this point, I couldn't figure out how to plot to a projectile motion under this forces. In other words, I'm trying to compare the projectile motion of the mass under three circumstances




  1. Mass only under the effect of gravity

  2. Mass under the effect of gravity and air resistance

  3. Mass under the effect of gravity, air resistance, and buoyancy






python numpy matplotlib physics







share|improve this question




















  • 1





    You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

    – Bazingaa
    Jan 2 at 8:19













  • From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

    – Mstaino
    Jan 2 at 13:05











  • I edited my question and tried to explain more clearly.

    – Birkan Emrem
    Jan 2 at 13:28














1












1








1








I am trying to make a plot of a projectile motion of a mass which is under the effect of gravitational, buoyancy and drag force. Basically, I want to show effects of the buoyancy and drag force on flight distance, flight time and velocity change on plotting.



import matplotlib.pyplot as plt

import numpy as np



V_initial = 30 # m/s

theta = np.pi/6 # 30

g = 3.711

m =1

C = 0.47

r = 0.5

S = np.pi*pow(r, 2)

ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)

t = np.linspace(0, t_flight, 200)



# Drag force

Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)



# Buoyant Force

Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))



x_loc = 

y_loc = 



for time in t:

    x = V_initial*time*np.cos(theta)

        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)

    x_loc.append(x)

    y_loc.append(y)



x_vel = 

y_vel = 

for time in t:

    vx = V_initial*np.cos(theta)

    vy = V_initial*np.sin(theta) - g*time

    x_vel.append(vx)

    y_vel.append(vy)





v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]



tau = 

for velocity in v_ch:

    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)

        tau.append(Ft)



buoy = 

for velocity in v_ch:

    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

    buoy.append(Fb)


after this point, I couldn't figure out how to plot to a projectile motion under this forces. In other words, I'm trying to compare the projectile motion of the mass under three circumstances




  1. Mass only under the effect of gravity

  2. Mass under the effect of gravity and air resistance

  3. Mass under the effect of gravity, air resistance, and buoyancy






python numpy matplotlib physics







share|improve this question
















I am trying to make a plot of a projectile motion of a mass which is under the effect of gravitational, buoyancy and drag force. Basically, I want to show effects of the buoyancy and drag force on flight distance, flight time and velocity change on plotting.



import matplotlib.pyplot as plt

import numpy as np



V_initial = 30 # m/s

theta = np.pi/6 # 30

g = 3.711

m =1

C = 0.47

r = 0.5

S = np.pi*pow(r, 2)

ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)

t = np.linspace(0, t_flight, 200)



# Drag force

Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)



# Buoyant Force

Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))



x_loc = 

y_loc = 



for time in t:

    x = V_initial*time*np.cos(theta)

        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)

    x_loc.append(x)

    y_loc.append(y)



x_vel = 

y_vel = 

for time in t:

    vx = V_initial*np.cos(theta)

    vy = V_initial*np.sin(theta) - g*time

    x_vel.append(vx)

    y_vel.append(vy)





v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]



tau = 

for velocity in v_ch:

    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)

        tau.append(Ft)



buoy = 

for velocity in v_ch:

    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

    buoy.append(Fb)


after this point, I couldn't figure out how to plot to a projectile motion under this forces. In other words, I'm trying to compare the projectile motion of the mass under three circumstances




  1. Mass only under the effect of gravity

  2. Mass under the effect of gravity and air resistance

  3. Mass under the effect of gravity, air resistance, and buoyancy






python numpy matplotlib physics




python numpy matplotlib physics






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited Jan 3 at 2:52









William Miller

1,426317




1,426317










asked Jan 2 at 2:06









Birkan EmremBirkan Emrem

216




216








  • 1





    You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

    – Bazingaa
    Jan 2 at 8:19













  • From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

    – Mstaino
    Jan 2 at 13:05











  • I edited my question and tried to explain more clearly.

    – Birkan Emrem
    Jan 2 at 13:28














  • 1





    You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

    – Bazingaa
    Jan 2 at 8:19













  • From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

    – Mstaino
    Jan 2 at 13:05











  • I edited my question and tried to explain more clearly.

    – Birkan Emrem
    Jan 2 at 13:28








1




1





You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

– Bazingaa
Jan 2 at 8:19







You gave a lot of undocumented code. You didn't specify in the question what variable you want to plot against which variable. You didn't tell us how the desired plot looks like (1d, 2d, 3d). You didn't tell us what exactly is the problem that you faced. How do you expect the readers to help you?

– Bazingaa
Jan 2 at 8:19















From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

– Mstaino
Jan 2 at 13:05





From what I understand your problem is twofold: solve (numerically) a differential equation, and plot it. Please elaborate a little bit on what you are trying to accomplish on both ends, and what you tried. You may also want to check scipy.integrate.odeint examples (we can help you if you specify the question better).

– Mstaino
Jan 2 at 13:05













I edited my question and tried to explain more clearly.

– Birkan Emrem
Jan 2 at 13:28





I edited my question and tried to explain more clearly.

– Birkan Emrem
Jan 2 at 13:28












2 Answers
2






active

oldest

votes


















2














You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.



F_x = tau_x 

F_y = tau_y + bouyancy + gravity


Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by 



v_x = v_x + (F_x / (2 * m)) * dt

v_y = v_y + (F_y / (2 * m)) * dt


So the x and y positions, r_x and r_y, at any time t are given by the summation of 



r_x = r_x + v_x * dt

r_y = r_y + v_y * dt


In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.



r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2

r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2


The entire calculation can actually be done in two lines, 



r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2

r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2


Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.



The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.



Edit



In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added



Complete example (edited) 



import numpy as np

import matplotlib.pyplot as plt



def projectile(V_initial, theta, bouyancy=True, drag=True):

    g = 9.81

    m = 1

    C = 0.47

    r = 0.5

    S = np.pi*pow(r, 2)

    ro_mars = 0.0175



    time = np.linspace(0, 100, 10000)

    tof = 0.0

    dt = time[1] - time[0]

    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))

    gravity = -g * m

    V_ix = V_initial * np.cos(theta)

    V_iy = V_initial * np.sin(theta)

    v_x = V_ix

    v_y = V_iy

    r_x = 0.0

    r_y = 0.0

    r_xs = list()

    r_ys = list()

    r_xs.append(r_x)

    r_ys.append(r_y)

    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.

    # Just make sure you use sufficiently small dt (dt is change in time between steps)

    for t in time:

        F_x = 0.0

        F_y = 0.0

        if (bouyancy == True):

            F_y = F_y + bouy

        if (drag == True):

            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)

            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)

        F_y = F_y + gravity



        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2

        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2

        v_x = v_x + (F_x / m) * dt

        v_y = v_y + (F_y / m) * dt

        if (r_y >= 0.0):

            r_xs.append(r_x)

            r_ys.append(r_y)

        else:

            tof = t

            r_xs.append(r_x)

            r_ys.append(r_y)

            break



    return r_xs, r_ys, tof



v = 30

theta = np.pi/4



fig = plt.figure(figsize=(8,4), dpi=300)

r_xs, r_ys, tof = projectile(v, theta, True, True)

plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")

r_xs, r_ys, tof = projectile(v, theta, False, True)

plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")

r_xs, r_ys, tof = projectile(v, theta, False, False)

plt.plot(r_xs, r_ys, 'k:', label="Gravity")

plt.title("Trajectory", fontsize=14)

plt.xlabel("Displacement in x-direction (m)")

plt.ylabel("Displacement in y-direction (m)")

plt.ylim(bottom=0.0)

plt.legend()

plt.show()




Note that this preserves and returns the time-of-flight in the variable tof.






share|improve this answer


























  • Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

    – Birkan Emrem
    Jan 15 at 7:54



















0














Using vector notation, and odeint.



import numpy as np

from scipy.integrate import odeint

import scipy.constants as SPC

import matplotlib.pyplot as plt



V_initial = 30 # m/s

theta = np.pi/6 # 30

g = 3.711

m = 1 # I assume this is your mass

C = 0.47

r = 0.5

ro_mars = 0.0175



t_flight = 2*(V_initial*np.sin(theta)/g)

t = np.linspace(0, t_flight, 200)



pos0 = [0, 0]

v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]



def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):

    x, y, x_prime, y_prime = vector



    # volume and surface

    V = np.pi * 4/3 * r**3

    S = np.pi*pow(r, 2)



    # net weight bouyancy

    if apply_bouyancy:

        Fb = (ro_mars * V - m) * g *np.array([0,1])

    else:

        Fb = -m  * g * np.array([0,1])



    # velocity vector

    v = np.array([x_prime, y_prime])



    # drag force - corrected to be updated based on current velocity

#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)

    if apply_resistance:

        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)

    else:

        Ft = np.array([0, 0])



    # resulting acceleration

    x_prime2, y_prime2 = (Fb + Ft) / m



    return x_prime, y_prime, x_prime2, y_prime2



sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))

plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')

plt.legend(loc='best')

plt.xlabel('x')

plt.ylabel('y')

plt.grid()

plt.show()


Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.



Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.



To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)






share|improve this answer























    Your Answer






    StackExchange.ifUsing("editor", function () {
    StackExchange.using("externalEditor", function () {
    StackExchange.using("snippets", function () {
    StackExchange.snippets.init();
    });
    });
    }, "code-snippets");

    StackExchange.ready(function() {
    var channelOptions = {
    tags: "".split(" "),
    id: "1"
    };
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function() {
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled) {
    StackExchange.using("snippets", function() {
    createEditor();
    });
    }
    else {
    createEditor();
    }
    });

    function createEditor() {
    StackExchange.prepareEditor({
    heartbeatType: 'answer',
    autoActivateHeartbeat: false,
    convertImagesToLinks: true,
    noModals: true,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: 10,
    bindNavPrevention: true,
    postfix: "",
    imageUploader: {
    brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
    contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
    allowUrls: true
    },
    onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    });


    }
    });














    draft saved

    draft discarded


















    StackExchange.ready(
    function () {
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f54000456%2fhow-to-plot-the-motion-of-a-projectile-under-the-effect-of-gravity-buoyancy-and%23new-answer', 'question_page');
    }
    );

    Post as a guest















    Required, but never shown

























    2 Answers
    2






    active

    oldest

    votes








    2 Answers
    2






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes









    2














    You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.



    F_x = tau_x 
    
F_y = tau_y + bouyancy + gravity


    Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by 



    v_x = v_x + (F_x / (2 * m)) * dt
    
v_y = v_y + (F_y / (2 * m)) * dt


    So the x and y positions, r_x and r_y, at any time t are given by the summation of 



    r_x = r_x + v_x * dt
    
r_y = r_y + v_y * dt


    In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.



    r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2


    The entire calculation can actually be done in two lines, 



    r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2


    Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.



    The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.



    Edit



    In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added



    Complete example (edited) 



    import numpy as np
    
import matplotlib.pyplot as plt
    

    
def projectile(V_initial, theta, bouyancy=True, drag=True):
    
    g = 9.81
    
    m = 1
    
    C = 0.47
    
    r = 0.5
    
    S = np.pi*pow(r, 2)
    
    ro_mars = 0.0175
    

    
    time = np.linspace(0, 100, 10000)
    
    tof = 0.0
    
    dt = time[1] - time[0]
    
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    
    gravity = -g * m
    
    V_ix = V_initial * np.cos(theta)
    
    V_iy = V_initial * np.sin(theta)
    
    v_x = V_ix
    
    v_y = V_iy
    
    r_x = 0.0
    
    r_y = 0.0
    
    r_xs = list()
    
    r_ys = list()
    
    r_xs.append(r_x)
    
    r_ys.append(r_y)
    
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    
    for t in time:
    
        F_x = 0.0
    
        F_y = 0.0
    
        if (bouyancy == True):
    
            F_y = F_y + bouy
    
        if (drag == True):
    
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
    
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
    
        F_y = F_y + gravity
    

    
        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
    
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
    
        v_x = v_x + (F_x / m) * dt
    
        v_y = v_y + (F_y / m) * dt
    
        if (r_y >= 0.0):
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
        else:
    
            tof = t
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
            break
    

    
    return r_xs, r_ys, tof
    

    
v = 30
    
theta = np.pi/4
    

    
fig = plt.figure(figsize=(8,4), dpi=300)
    
r_xs, r_ys, tof = projectile(v, theta, True, True)
    
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, True)
    
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, False)
    
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
    
plt.title("Trajectory", fontsize=14)
    
plt.xlabel("Displacement in x-direction (m)")
    
plt.ylabel("Displacement in y-direction (m)")
    
plt.ylim(bottom=0.0)
    
plt.legend()
    
plt.show()




    Note that this preserves and returns the time-of-flight in the variable tof.






    share|improve this answer


























    • Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

      – Birkan Emrem
      Jan 15 at 7:54
















    2














    You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.



    F_x = tau_x 
    
F_y = tau_y + bouyancy + gravity


    Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by 



    v_x = v_x + (F_x / (2 * m)) * dt
    
v_y = v_y + (F_y / (2 * m)) * dt


    So the x and y positions, r_x and r_y, at any time t are given by the summation of 



    r_x = r_x + v_x * dt
    
r_y = r_y + v_y * dt


    In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.



    r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2


    The entire calculation can actually be done in two lines, 



    r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2


    Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.



    The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.



    Edit



    In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added



    Complete example (edited) 



    import numpy as np
    
import matplotlib.pyplot as plt
    

    
def projectile(V_initial, theta, bouyancy=True, drag=True):
    
    g = 9.81
    
    m = 1
    
    C = 0.47
    
    r = 0.5
    
    S = np.pi*pow(r, 2)
    
    ro_mars = 0.0175
    

    
    time = np.linspace(0, 100, 10000)
    
    tof = 0.0
    
    dt = time[1] - time[0]
    
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    
    gravity = -g * m
    
    V_ix = V_initial * np.cos(theta)
    
    V_iy = V_initial * np.sin(theta)
    
    v_x = V_ix
    
    v_y = V_iy
    
    r_x = 0.0
    
    r_y = 0.0
    
    r_xs = list()
    
    r_ys = list()
    
    r_xs.append(r_x)
    
    r_ys.append(r_y)
    
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    
    for t in time:
    
        F_x = 0.0
    
        F_y = 0.0
    
        if (bouyancy == True):
    
            F_y = F_y + bouy
    
        if (drag == True):
    
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
    
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
    
        F_y = F_y + gravity
    

    
        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
    
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
    
        v_x = v_x + (F_x / m) * dt
    
        v_y = v_y + (F_y / m) * dt
    
        if (r_y >= 0.0):
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
        else:
    
            tof = t
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
            break
    

    
    return r_xs, r_ys, tof
    

    
v = 30
    
theta = np.pi/4
    

    
fig = plt.figure(figsize=(8,4), dpi=300)
    
r_xs, r_ys, tof = projectile(v, theta, True, True)
    
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, True)
    
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, False)
    
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
    
plt.title("Trajectory", fontsize=14)
    
plt.xlabel("Displacement in x-direction (m)")
    
plt.ylabel("Displacement in y-direction (m)")
    
plt.ylim(bottom=0.0)
    
plt.legend()
    
plt.show()




    Note that this preserves and returns the time-of-flight in the variable tof.






    share|improve this answer


























    • Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

      – Birkan Emrem
      Jan 15 at 7:54














    2












    2








    2







    You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.



    F_x = tau_x 
    
F_y = tau_y + bouyancy + gravity


    Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by 



    v_x = v_x + (F_x / (2 * m)) * dt
    
v_y = v_y + (F_y / (2 * m)) * dt


    So the x and y positions, r_x and r_y, at any time t are given by the summation of 



    r_x = r_x + v_x * dt
    
r_y = r_y + v_y * dt


    In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.



    r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2


    The entire calculation can actually be done in two lines, 



    r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2


    Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.



    The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.



    Edit



    In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added



    Complete example (edited) 



    import numpy as np
    
import matplotlib.pyplot as plt
    

    
def projectile(V_initial, theta, bouyancy=True, drag=True):
    
    g = 9.81
    
    m = 1
    
    C = 0.47
    
    r = 0.5
    
    S = np.pi*pow(r, 2)
    
    ro_mars = 0.0175
    

    
    time = np.linspace(0, 100, 10000)
    
    tof = 0.0
    
    dt = time[1] - time[0]
    
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    
    gravity = -g * m
    
    V_ix = V_initial * np.cos(theta)
    
    V_iy = V_initial * np.sin(theta)
    
    v_x = V_ix
    
    v_y = V_iy
    
    r_x = 0.0
    
    r_y = 0.0
    
    r_xs = list()
    
    r_ys = list()
    
    r_xs.append(r_x)
    
    r_ys.append(r_y)
    
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    
    for t in time:
    
        F_x = 0.0
    
        F_y = 0.0
    
        if (bouyancy == True):
    
            F_y = F_y + bouy
    
        if (drag == True):
    
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
    
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
    
        F_y = F_y + gravity
    

    
        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
    
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
    
        v_x = v_x + (F_x / m) * dt
    
        v_y = v_y + (F_y / m) * dt
    
        if (r_y >= 0.0):
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
        else:
    
            tof = t
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
            break
    

    
    return r_xs, r_ys, tof
    

    
v = 30
    
theta = np.pi/4
    

    
fig = plt.figure(figsize=(8,4), dpi=300)
    
r_xs, r_ys, tof = projectile(v, theta, True, True)
    
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, True)
    
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, False)
    
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
    
plt.title("Trajectory", fontsize=14)
    
plt.xlabel("Displacement in x-direction (m)")
    
plt.ylabel("Displacement in y-direction (m)")
    
plt.ylim(bottom=0.0)
    
plt.legend()
    
plt.show()




    Note that this preserves and returns the time-of-flight in the variable tof.






    share|improve this answer















    You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.



    F_x = tau_x 
    
F_y = tau_y + bouyancy + gravity


    Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by 



    v_x = v_x + (F_x / (2 * m)) * dt
    
v_y = v_y + (F_y / (2 * m)) * dt


    So the x and y positions, r_x and r_y, at any time t are given by the summation of 



    r_x = r_x + v_x * dt
    
r_y = r_y + v_y * dt


    In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.



    r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2


    The entire calculation can actually be done in two lines, 



    r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
    
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2


    Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.



    The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.



    Edit



    In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added



    Complete example (edited) 



    import numpy as np
    
import matplotlib.pyplot as plt
    

    
def projectile(V_initial, theta, bouyancy=True, drag=True):
    
    g = 9.81
    
    m = 1
    
    C = 0.47
    
    r = 0.5
    
    S = np.pi*pow(r, 2)
    
    ro_mars = 0.0175
    

    
    time = np.linspace(0, 100, 10000)
    
    tof = 0.0
    
    dt = time[1] - time[0]
    
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    
    gravity = -g * m
    
    V_ix = V_initial * np.cos(theta)
    
    V_iy = V_initial * np.sin(theta)
    
    v_x = V_ix
    
    v_y = V_iy
    
    r_x = 0.0
    
    r_y = 0.0
    
    r_xs = list()
    
    r_ys = list()
    
    r_xs.append(r_x)
    
    r_ys.append(r_y)
    
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    
    for t in time:
    
        F_x = 0.0
    
        F_y = 0.0
    
        if (bouyancy == True):
    
            F_y = F_y + bouy
    
        if (drag == True):
    
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
    
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
    
        F_y = F_y + gravity
    

    
        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
    
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
    
        v_x = v_x + (F_x / m) * dt
    
        v_y = v_y + (F_y / m) * dt
    
        if (r_y >= 0.0):
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
        else:
    
            tof = t
    
            r_xs.append(r_x)
    
            r_ys.append(r_y)
    
            break
    

    
    return r_xs, r_ys, tof
    

    
v = 30
    
theta = np.pi/4
    

    
fig = plt.figure(figsize=(8,4), dpi=300)
    
r_xs, r_ys, tof = projectile(v, theta, True, True)
    
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, True)
    
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
    
r_xs, r_ys, tof = projectile(v, theta, False, False)
    
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
    
plt.title("Trajectory", fontsize=14)
    
plt.xlabel("Displacement in x-direction (m)")
    
plt.ylabel("Displacement in y-direction (m)")
    
plt.ylim(bottom=0.0)
    
plt.legend()
    
plt.show()




    Note that this preserves and returns the time-of-flight in the variable tof.







    share|improve this answer














    share|improve this answer



    share|improve this answer








    edited Feb 14 at 3:58

























    answered Jan 2 at 8:40









    William MillerWilliam Miller

    1,426317




    1,426317













    • Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

      – Birkan Emrem
      Jan 15 at 7:54



















    • Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

      – Birkan Emrem
      Jan 15 at 7:54

















    Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

    – Birkan Emrem
    Jan 15 at 7:54





    Hi William, what is the value of the gravity parameter in function it seems it's not defined and because of that it does not give the graph you found.

    – Birkan Emrem
    Jan 15 at 7:54













    0














    Using vector notation, and odeint.



    import numpy as np
    
from scipy.integrate import odeint
    
import scipy.constants as SPC
    
import matplotlib.pyplot as plt
    

    
V_initial = 30 # m/s
    
theta = np.pi/6 # 30
    
g = 3.711
    
m = 1 # I assume this is your mass
    
C = 0.47
    
r = 0.5
    
ro_mars = 0.0175
    

    
t_flight = 2*(V_initial*np.sin(theta)/g)
    
t = np.linspace(0, t_flight, 200)
    

    
pos0 = [0, 0]
    
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]
    

    
def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
    
    x, y, x_prime, y_prime = vector
    

    
    # volume and surface
    
    V = np.pi * 4/3 * r**3
    
    S = np.pi*pow(r, 2)
    

    
    # net weight bouyancy
    
    if apply_bouyancy:
    
        Fb = (ro_mars * V - m) * g *np.array([0,1])
    
    else:
    
        Fb = -m  * g * np.array([0,1])
    

    
    # velocity vector
    
    v = np.array([x_prime, y_prime])
    

    
    # drag force - corrected to be updated based on current velocity
    
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
    
    if apply_resistance:
    
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
    
    else:
    
        Ft = np.array([0, 0])
    

    
    # resulting acceleration
    
    x_prime2, y_prime2 = (Fb + Ft) / m
    

    
    return x_prime, y_prime, x_prime2, y_prime2
    

    
sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
    
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
    
plt.legend(loc='best')
    
plt.xlabel('x')
    
plt.ylabel('y')
    
plt.grid()
    
plt.show()


    Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.



    Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.



    To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)






    share|improve this answer




























      0














      Using vector notation, and odeint.



      import numpy as np
      
from scipy.integrate import odeint
      
import scipy.constants as SPC
      
import matplotlib.pyplot as plt
      

      
V_initial = 30 # m/s
      
theta = np.pi/6 # 30
      
g = 3.711
      
m = 1 # I assume this is your mass
      
C = 0.47
      
r = 0.5
      
ro_mars = 0.0175
      

      
t_flight = 2*(V_initial*np.sin(theta)/g)
      
t = np.linspace(0, t_flight, 200)
      

      
pos0 = [0, 0]
      
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]
      

      
def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
      
    x, y, x_prime, y_prime = vector
      

      
    # volume and surface
      
    V = np.pi * 4/3 * r**3
      
    S = np.pi*pow(r, 2)
      

      
    # net weight bouyancy
      
    if apply_bouyancy:
      
        Fb = (ro_mars * V - m) * g *np.array([0,1])
      
    else:
      
        Fb = -m  * g * np.array([0,1])
      

      
    # velocity vector
      
    v = np.array([x_prime, y_prime])
      

      
    # drag force - corrected to be updated based on current velocity
      
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
      
    if apply_resistance:
      
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
      
    else:
      
        Ft = np.array([0, 0])
      

      
    # resulting acceleration
      
    x_prime2, y_prime2 = (Fb + Ft) / m
      

      
    return x_prime, y_prime, x_prime2, y_prime2
      

      
sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
      
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
      
plt.legend(loc='best')
      
plt.xlabel('x')
      
plt.ylabel('y')
      
plt.grid()
      
plt.show()


      Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.



      Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.



      To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)






      share|improve this answer


























        0












        0








        0







        Using vector notation, and odeint.



        import numpy as np
        
from scipy.integrate import odeint
        
import scipy.constants as SPC
        
import matplotlib.pyplot as plt
        

        
V_initial = 30 # m/s
        
theta = np.pi/6 # 30
        
g = 3.711
        
m = 1 # I assume this is your mass
        
C = 0.47
        
r = 0.5
        
ro_mars = 0.0175
        

        
t_flight = 2*(V_initial*np.sin(theta)/g)
        
t = np.linspace(0, t_flight, 200)
        

        
pos0 = [0, 0]
        
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]
        

        
def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
        
    x, y, x_prime, y_prime = vector
        

        
    # volume and surface
        
    V = np.pi * 4/3 * r**3
        
    S = np.pi*pow(r, 2)
        

        
    # net weight bouyancy
        
    if apply_bouyancy:
        
        Fb = (ro_mars * V - m) * g *np.array([0,1])
        
    else:
        
        Fb = -m  * g * np.array([0,1])
        

        
    # velocity vector
        
    v = np.array([x_prime, y_prime])
        

        
    # drag force - corrected to be updated based on current velocity
        
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
        
    if apply_resistance:
        
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
        
    else:
        
        Ft = np.array([0, 0])
        

        
    # resulting acceleration
        
    x_prime2, y_prime2 = (Fb + Ft) / m
        

        
    return x_prime, y_prime, x_prime2, y_prime2
        

        
sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
        
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
        
plt.legend(loc='best')
        
plt.xlabel('x')
        
plt.ylabel('y')
        
plt.grid()
        
plt.show()


        Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.



        Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.



        To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)






        share|improve this answer













        Using vector notation, and odeint.



        import numpy as np
        
from scipy.integrate import odeint
        
import scipy.constants as SPC
        
import matplotlib.pyplot as plt
        

        
V_initial = 30 # m/s
        
theta = np.pi/6 # 30
        
g = 3.711
        
m = 1 # I assume this is your mass
        
C = 0.47
        
r = 0.5
        
ro_mars = 0.0175
        

        
t_flight = 2*(V_initial*np.sin(theta)/g)
        
t = np.linspace(0, t_flight, 200)
        

        
pos0 = [0, 0]
        
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]
        

        
def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
        
    x, y, x_prime, y_prime = vector
        

        
    # volume and surface
        
    V = np.pi * 4/3 * r**3
        
    S = np.pi*pow(r, 2)
        

        
    # net weight bouyancy
        
    if apply_bouyancy:
        
        Fb = (ro_mars * V - m) * g *np.array([0,1])
        
    else:
        
        Fb = -m  * g * np.array([0,1])
        

        
    # velocity vector
        
    v = np.array([x_prime, y_prime])
        

        
    # drag force - corrected to be updated based on current velocity
        
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
        
    if apply_resistance:
        
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
        
    else:
        
        Ft = np.array([0, 0])
        

        
    # resulting acceleration
        
    x_prime2, y_prime2 = (Fb + Ft) / m
        

        
    return x_prime, y_prime, x_prime2, y_prime2
        

        
sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
        
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
        
plt.legend(loc='best')
        
plt.xlabel('x')
        
plt.ylabel('y')
        
plt.grid()
        
plt.show()


        Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.



        Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.



        To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)







        share|improve this answer












        share|improve this answer



        share|improve this answer










        answered Jan 2 at 18:25









        MstainoMstaino

        2,0221413




        2,0221413






























            draft saved

            draft discarded




















































            Thanks for contributing an answer to Stack Overflow!


            • Please be sure to answer the question. Provide details and share your research!

            But avoid …



            • Asking for help, clarification, or responding to other answers.

            • Making statements based on opinion; back them up with references or personal experience.


            To learn more, see our tips on writing great answers.




            draft saved


            draft discarded














            StackExchange.ready(
            function () {
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f54000456%2fhow-to-plot-the-motion-of-a-projectile-under-the-effect-of-gravity-buoyancy-and%23new-answer', 'question_page');
            }
            );

            Post as a guest















            Required, but never shown





















































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown

































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown







            Popular posts from this blog

            MongoDB - Not Authorized To Execute Command

            Image
            .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty{ height:90px;width:728px;box-sizing:border-box; } 0 I have successfully enabled authorization on MongoDB and I have created an account on the admin database and then I created an account for my database called test. The following connection string to connect to my test database works successfully: mongo --host 192.168.17.52 --port 27017 -u user1 -p password --authenticationDatabase test Only problem I have now is, I cannot execute commands such as: show dbs. I get the following error when I try to do so: "errmsg" : "not authorized on admin to execute command { listDatabases: 1.0, lsid: { id: UUID("a1d5bc0d-bc58-485e-b232-270758a89455") }, $db: "admin...
            Read more

            How to fix TextFormField cause rebuild widget in Flutter

            Image
            .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty{ height:90px;width:728px;box-sizing:border-box; } 0 I've got the same issue at #9280. I tried all solutions from that topic with/without GlobalKey for Scaffold and Form but the problem wasn't fixed. I tried to build a authentication form with login and signup. But the keyboard doesn't show when I tap the TextFormField at the first time. Then I try to input some letter from physical keyboard, the widget is rebuilt. Like this: login screen Nothing to be log when I run flutter run --verbose, so I logged manually. Every tapping on field, the widget called dispose -> init -> build. [√] Flutter (Channel beta, v1.0.0, on Microsoft Windows [Version 10.0.17763.195], ...
            Read more

            in spring boot 2.1 many test slices are not allowed anymore due to multiple @BootstrapWith

            Image
            1 I tried to upgrade a yummy sandwich made of two test slices (@JsonTest and @JdbcTest in my case, crunchy test code in between) adding spring boot 2.1 flavour to it. But it seems it was not much of a success. I cannot annotate my tests with many @...Test since they are now each bringing their own XxxTestContextBootstrapper. It used to work when they all used same SpringBootTestContextBootstrapper. @RunWith(SpringRunner.class) @JdbcTest @JsonTest public class Test { @Test public void test() { System.out.printn("Hello, World !"); } } The error I get from BootstrapUtils is illegalStateException : Configuration error: found multiple declarations of @BootstrapWith for test class I understand I might be doing something wrong here but is there an easy way I could load both Json and Jdbc contex...
            Read more