0

I am trying to model firing a projectile from a slingshot.

This is my code:

from pylab import * 
import numpy as np
from scipy.integrate import odeint
import seaborn
## set initial conditions and parameters
g = 9.81            # acceleration due to gravity
th = 30            # set launch angle
th = th * np.pi/180.   # convert launch angle to radians
v0 = 10.0           # set initial speed
c = 0.5             # controls strength of air drag
d = 0.02 # diameter of the spherical rock
A = pi * (d/2)**2 #the cross-sectional area of the spherical rock
ro = 1.2041 #the density of the medium we are perfoming the launch in
m = 0.01 #mass

x0=0                # specify initial conditions
y0=0
vx0 = v0*sin(th)
vy0 = v0*cos(th)

## defining our model
def slingshot_model(state,time):
    z = zeros(4)    # create array to hold z vector
    z[0] = state[2] # z[0] = x component of velocity
    z[1] = state[3] # z[1] = y component of velocity
    z[2] = - c*A*ro/2*m*sqrt(z[0]**2 + z[1]**2)*z[0]         # z[2] = acceleration in x direction
    z[3] = -g/m - c*A*ro/2*m*sqrt(z[0]**2 + z[1]**2)*z[1]       # z[3] = acceleration in y direction
    return z

## set initial state vector and time array
X0 = [x0, y0, vx0, vy0]        # set initial state of the system
t0 = 0
tf = 4 #final time
tau = 0.05 #time step


# create time array starting at t0, ending at tf with a spacing tau
t = arange(t0,tf,tau)   

## solve ODE using odeint
X = odeint(slingshot_model,X0,t) # returns an 2-dimensional array with the 
                        # first index specifying the time and the
                        # second index specifying the component of
                        # the state vector

# putting ':' as an index specifies all of the elements for
# that index so x, y, vx, and vy are arrays at times specified 
# in the time array
x = X[:,0]  
y = X[:,1] 
vx = X[:,2] 
vy = X[:,3]

plt.rcParams['figure.figsize'] = [10, 10]
plot(x,y)

But it gives me this plot that doesn't make sense to me:

enter image description here

What am I missing? The values shouldn't come out like they do, but for the life of me I can't see why.

It is probably something trivial, but I have been staring at this too long, so I figured bringing in a fresh set of eyes is the best course of action.

nedimk994
  • 9
  • 2
  • It is a wonderful parabola trajectory with a seemingly horizontal launch angle. What it wrong with it? – Reblochon Masque May 29 '18 at 05:30
  • Hello, thanks for the reply, As part of my calculations I also have to do the following: You shoot from a height of 0m. What is the range of your sling-shot on an even plane (no obstacles, elevation etc...) Find out when and where the rock falls down, that is, what are t and x when y=0. Play a little with the angle and try to maximize the range! But how can I find t and x when y goes from 0 to -8000 on the plot. This is what bugs me. – nedimk994 May 29 '18 at 05:36
  • You need to shoot up a bit. – Reblochon Masque May 29 '18 at 05:38
  • @ReblochonMasque He is shooting up at an angle of 60 deg wrt to X axis – AGN Gazer May 29 '18 at 05:39
  • `th = 30 # set launch angle`... – Reblochon Masque May 29 '18 at 05:42
  • @ReblochonMasque 30 with regard to what? – AGN Gazer May 29 '18 at 05:43
  • Launch angle... wrt the imperial battlecruiser azimuth? Poor choices of coordinate systems and/or confusion of sin with cos and vice versa have a cost! – Reblochon Masque May 29 '18 at 05:44
  • Yes, that is my point: launch_angle of projectile in an x, y coordinate system, is assumed to be wrt x, unless otherwise specified. – Reblochon Masque May 29 '18 at 05:49
  • The point is: if the angle is defined "as usual" with regard to the positive direction of the X axis then OP computations for the projection of the velocity are wrong. – AGN Gazer May 29 '18 at 05:49
  • Yes, they are, and you pointed this out in your answer. – Reblochon Masque May 29 '18 at 05:50
  • All I was saying is that his (incorrect) computations are equivalent to a launch at 60 deg using correct computations. – AGN Gazer May 29 '18 at 05:52
  • You first should get correct picture for a simpler problem - no drag. Therefore set `c = 0` and check what you get. And then you will realize that `z[3] = -g/m` (which you define as `z[3] = acceleration in y direction`) **makes no sense at all**: your acceleration is proportional to `g` (already an acceleration) divided by mass! – AGN Gazer May 29 '18 at 05:55

1 Answers1

1

I think there are at least two major problems with your computations:

  1. Usually angle is defined with regard to the X-axis. Therefore

    vx0 = v0*cos(th) # not sin
vy0 = v0*sin(th) # not cos

  2. Most importantly, why are you dividing acceleration of the free fall g by the mass? (see z[3] = -g/m...) This makes no sense to me. DO NOT divide by mass!

EDIT:

  1. Based on your comment and linked formulae, it is clear that your code also suffers from a third mistake: air drag terms should be inverse-proportional to mass:

enter image description here

AGN Gazer
  • 8,025
  • 2
  • 27
  • 45
  • Thank you for your answer, I have edited my equations as per your instructions and now the plots make sense. As for dividing by mass these are the equations I was given y my professor: https://i.stack.imgur.com/yueH6.png – nedimk994 May 29 '18 at 05:54
  • @nedimk994 Then you have another error in your code. See my last edit. – AGN Gazer May 29 '18 at 06:01
  • @nedimk994 Please consider accepting this answer if it solved your problem. – AGN Gazer May 30 '18 at 02:30