Physics with Python. Full article at http://www.firsttimeprogrammer.blogspot.com/2014/07/physics-with-python.html

physics_py.py

import math
import numpy as np
import matplotlib.pyplot as plt
import pylab
import os

# this is optional, if you want to save pictures of the trajectory
# later, then choose a folder where to save images
os.chdir("path")

# Acceleration of gravity = 9.80665 m/s**2

class Projectile_motion(object):
# s0 is horizontal starting space, we assume that the projectile starts at y = 0
# v0 is the initial velocity
# theta is the launch angle in radiants
# v0x is the horizontal component of velocity
# v0y is the vertical component of velocity
    def __init__(self,s0,v0,theta):
        self.s0 = s0
        self.v0 = v0
        self.theta = theta
        self.v0x = v0*math.cos(theta)
        self.v0y = v0*math.sin(theta)

    # The representation of an instance of this class
    def __repr__(self):
        return ("The projectile starts at %s with an initial velocity of %s meters per second") %(round(self.s0,2),round(self.v0,2))
    
    # The maximum range of the instance
    def range_(self):
        space = (2* (self.v0)**2 * math.sin(self.theta) * math.cos(self.theta))/(9.81)
        space_km = space/1000
        return "The range is %s m, that's to say: %s km." %(space,space_km)
    
    # The maximum achievable height
    def max_height(self):
        height = ((self.v0)**2 * math.sin(self.theta)**2)/(a*2)
        height_km = height/1000
        return "Maximum achievable height is %s m, that's to say: %s km." %(height,height_km)
    
    # For how long does the projectile stay up in the air?
    def time_in_air(self):
        time = (2*self.v0*math.sin(self.theta))/a
        return time
    
    # This function gives the speed at each time t
    def speed(self,t):
        if t > self.time_in_air():
            return "Projectile has already landed, please reduce t. Projectile landed at %s" %(self.time_in_air())
        y = self.v0y - a*t
        x = self.v0x
        vel = math.sqrt(x**2 + y**2)
        if y > 0:
            return "At %s seconds: horizontal speed (constant): %s m/s. Vertical speed %sm/s upwards. Combined speed: %sm/s" %(round(t,0),round(x,0),round(y,2),round(vel,2))
        elif y == 0:
            return "At %s seconds: horizontal speed (constant): %s m/s. Vertical speed 0. Combined speed: %sm/s" %(round(t,0),round(x,0),round(y,2),round(vel,2))
        else:
            return "At %s seconds: horizontal speed (constant): %s m/s. Vertical speed %sm/s downwards. Combined speed: %sm/s" %(round(t,0),round(x,0),round(y,2),round(vel,2))
    
    # How far can it travel?
    def how_far(self,t):
        if t > self.time_in_air():
            return "Projectile has already landed, please reduce t. Projectile landed at %s" %(self.time_in_air())
        x = self.v0x*t - self.s0
        x_km = x/1000
        return "Horizontal space travelled in %s seconds is equal to %s m or, %s km" %(t,x,x_km)
    
    # This function returns the position at every time t
    def position(self,t):
        if t > self.time_in_air():
            return "Projectile has already landed, please reduce t. Projectile landed at %s" %(self.time_in_air())
        x = self.s0 + self.v0x*t
        y = self.s0 + self.v0y*t - 0.5*a*t**2
        position_t = []
        position_t.append(x)
        position_t.append(y)
        return position_t
    
    # This is the most interesting function, it returns the graph of the trajectory
    # Optionally you can save it or save a copy of each stage of graphing
    def show_trajectory(self):
        i = 0
        x = []
        y = []
        while i <= self.time_in_air():
            x.append(self.position(i)[0])
            y.append(self.position(i)[1])
            i += 1
            plot = plt.plot(x,y,"ro")
            # If you decide to print out a copy of the graph at each stage then remove
            # "#"from the next line of code. This will save for each i a .png picture in the path specified above with os.chdir()
            #pylab.savefig("hey"+"_%s_"%i+".png")
            plt.xlabel("X")
            plt.ylabel("Y")
            plt.show()p = Projectile_motion(0,600,np.pi/3)