Demo entry 1967976

trajectory calculator

   

Submitted by anonymous on Jun 20, 2015 at 11:41
Language: Python. Code size: 2.6 kB.

# -*- coding: utf-8 -*-
from numpy import *
from pylab import *

radius = 1.5e-3 #m
densityWater = 1000.0 #kg/m3
mass = 4.0/3.0*pi*radius**3*densityWater #kg
Const = 0.6
densityAir = 1.2 #kg/m3
Area = pi*radius**2
print radius,mass,Area

theta=pi/3. #rad
vel=1.0 #m/s

nt=51
timeTot=1.0 #s
dt=timeTot/(nt-1)
time=arange(nt)*dt

gx=0.0
gy=-10.0 #m/s2

aGx = ones(nt-2) * gx; aGy = ones(nt-2) * gy
aDx = zeros(nt-2); aDy = zeros(nt-2)
vx = zeros(nt-1); vy = zeros(nt-1)
vx[0]=vel*cos(theta)
vy[0]=vel*sin(theta)

for it in arange(nt-2):
    vMag2=vx[it]**2+vy[it]**2
    accDrag=0.5*Const*densityAir*Area*vMag2/mass
    if accDrag>0.:
        aDx[it]=-accDrag*vx[it]/vMag2**0.5
        aDy[it]=-accDrag*vy[it]/vMag2**0.5
    else:
        aDx[it]=0
        aDy[it]=0
    vx[it+1] = vx[it]+(aGx[it]+aDx[it])*dt
    vy[it+1] = vy[it]+(aGy[it]+aDy[it])*dt
    print vMag2,accDrag,aDx[it],aDy[it],vx[it],vx[it+1],vy[it],vy[it+1]

pxApprox = zeros(nt)
pyApprox = zeros(nt)
for it in arange(nt-1):
    pxApprox[it+1] = pxApprox[it]+vx[it]*dt
    pyApprox[it+1] = pyApprox[it]+vy[it]*dt
                
px = zeros(nt)
py = zeros(nt)
for it in arange(nt-1):
    if it<nt-2:
        px[it+1] = px[it] + 0.5*(vx[it]+vx[it+1])*dt
        py[it+1] = py[it] + 0.5*(vy[it]+vy[it+1])*dt
    else:
        px[it+1] = px[it]+vx[it]*dt
        py[it+1] = py[it]+vy[it]*dt

#benchmark with the analytical solution    
pxBench = zeros(nt)
pyBench = zeros(nt)
for it in arange(nt):
    pxBench[it]=vel*cos(theta)*time[it]
    pyBench[it]=vel*sin(theta)*time[it]-0.5*10.0*time[it]**2

#Kinetic Energy
KE = zeros(nt-1)
WorkDrag = zeros(nt-2)
for it in arange(nt-1):
    KE[it]=0.5*mass*(vx[it]**2+vy[it]**2)

#Gravitational Potential Energy
PE = zeros(nt)
for it in arange(nt):    
    PE[it]=mass*(-gy)*py[it]
    
#Mechanical Energy
ME = zeros(nt-1)
for it in arange(nt-1):
    ME[it]=KE[it]+0.5*(PE[it]+PE[it+1])    

#Work of the Drag Force
WDrag = zeros(nt-1)
for it in arange(nt-2):
    WDrag[it+1]=WDrag[it]-mass*(aDx[it]*(px[it+1]-px[it])+aDy[it]*(py[it+1]-py[it]))

#Total Energy
TotalEnergy = zeros(nt-1)
for it in arange(nt-1):
    TotalEnergy[it]=ME[it]+WDrag[it]

# plot everything
#l, = plot(pxApprox, pyApprox, 'g', label='Approximate')
#l, = plot(px, py, 'b', label='py vs px')
#l, = plot(pxBench, pyBench, 'r', label='Benchmark')

plot(time[0:nt-1],KE,label='KE')
plot(time,PE,label='PE')
plot(time[0:nt-1],ME,label='ME')
plot(time[0:nt-1],WDrag,label='WDrag')
plot(time[0:nt-1],TotalEnergy,label='TE')

legend(loc='lower left')
show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).