# Demo entry 2767809

Project 1

Submitted by anonymous on Sep 28, 2015 at 22:19
Language: Python. Code size: 5.4 kB.

```# -*- coding: utf-8 -*-
"""
Created on Sun Sep 27 15:37:01 2015

@author: RhyzP
"""

import numpy as np
import math
import pylab

#Finds the initial volume given initial height and s
def V_0(s, h):
r=h/s
return (1/3)*math.pi*h*r**2

def Euler(n,b,s,y_top=6.0,y_init=1.0):
#Use the user-specified b value to find y_out
y_out=y_top*b

#Initial height of the water
y=1.0

#Define the initial time and step size
t=0.5
step=0.5
#Find the initial radius using s and the initial height
r=y_init/s

#Use previous function to find intial volume
V0=V_0(s,y_init)

#Define the volume array which starts out empty and append initial volume
V_a=[]
V_a.append(V0)

#Define number of steps based on max value and step size
num=(n/step)
num=int(num)
#iterate through all the steps
for i in range(num):
#Have an if statement for whether the water is below y_out
if y>y_out:
#Euler's formula
Vt=V0+step*(3.5*(math.sin(t))**2-3.2*((3*V0)/(math.pi*r**2)-y_out))
r=(3*Vt/(math.pi*s))**(1.0/3.0)#The new radius of the volume of the liquid
y=r*s#The new height of the liquid
V_a.append(Vt)#append new volume to volume array
V0=Vt
t+=step
else:
Vt=V0+step*(3.5*(math.sin(t))**2)
r=(3*Vt/(math.pi*s))**(1.0/3.0)
y=r*s
V_a.append(Vt)
V0=Vt
t+=step
return V_a

#Makes a graph for five sets of data
def grapher(x, y, xlabel, ylabel, title):
pylab.figure(figsize=(8,6), dpi=150)
pylab.subplot(1,1,1)
types="r-", "g-", "b-", "y-", "k-"
m=len(y)
for i in range(m):
pylab.plot(x, y[i], types[i])
pylab.title(title)
pylab.xlabel(xlabel)
pylab.ylabel(ylabel)
pylab.grid(which='major',axis='both')

#returns an array of y values over time
def y_a(V_a,s):
n=0
for a in V_a:
a=((3.0*a*s**2)/(math.pi))**(1.0/3.0)
V_a[n]=a
n+=1
y_a=V_a
return y_a

#return the value of the flow at particular y values
def q_a(y_a,b, y_top=6.0):
q_a=[]
for a in y_a:
#a<y_out which changes based on b
if a<(6*b):
q_a.append(0)
else:
q_a.append(3.2*(a-b*6)**1.4)
return q_a

#Finds t_res for all cases
def t_res(q_a,V_a):
t_res=[]
n=0
for a in q_a:
if a==0:
t_res.append(V_a[n]/0.01)
n+=1
else:
t_res.append(V_a[n]/a)
n+=1
return t_res

def time_maker(initial, final, step):
fin=final+step
time=np.arange(initial, fin, step)
return time

# -*- coding: utf-8 -*-
"""
Created on Sun Sep 27 20:01:01 2015

@author: RhyzP
"""

import numpy as np
import math
import pylab
import Project1_modules as p1

""" Part A"""
def Part_A():
t_a=p1.time_maker(0,70,0.5)
b=0.25
s=1.6
V_a=p1.Euler(70,b,s)
y_a=p1.y_a(V_a,s)
q_a=p1.q_a(y_a,b)
t_res=p1.t_res(q_a,V_a)

pylab.figure(figsize=(8,6), dpi=80)
pylab.subplot(1,1,1)
pylab.plot(t_a, y_a, "ro:")
pylab.title('V(t) vs t')
pylab.xlabel('time')
pylab.ylabel('Volume')
pylab.grid(which='major',axis='both')

pylab.figure(figsize=(8,6), dpi=80)
pylab.subplot(1,1,1)
pylab.plot(t_a, t_res, "ro:")
pylab.title('V(t) vs t')
pylab.xlabel('Time')
pylab.ylabel('Residence Time')
pylab.grid(which='major',axis='both')

#Part_A()

"""Part B"""

def Part_B():
t_a=p1.time_maker(0,70,0.5)
s=1.6
y_out=np.arange(2,4.5,.5)
b=y_out/6.0
Y_a=[]
Q_a=[]
T_res=[]
for element in b:
V_a=p1.Euler(70,element,s)
y_a=p1.y_a(V_a,s)
Y_a.append(y_a)
q_a=p1.q_a(y_a,element)
Q_a.append(q_a)
t_res=p1.t_res(q_a,V_a)
T_res.append(t_res)

p1.grapher(t_a,Y_a,'Time','Height','H(t) vs t')

p1.grapher(t_a,Q_a,'Time','Outflow','Q(t) vs t')

p1.grapher(t_a,T_res,'Time','Residence Time','t_res(t) vs t')

Part_B()

"""Part C"""

def Part_C1():
t_a=p1.time_maker(0,70,0.5)
b=0.25
r_top=np.arange(3,5.5,.5)
s=6.0/r_top
Y_a=[]
Q_a=[]
T_res=[]
for element in s:
V_a=p1.Euler(70,b,element)
y_a=p1.y_a(V_a,element)
Y_a.append(y_a)
q_a=p1.q_a(y_a,b)
Q_a.append(q_a)
t_res=p1.t_res(q_a,V_a)
T_res.append(t_res)
p1.grapher(t_a,Y_a,'Time','Height','H(t) vs t')

p1.grapher(t_a,Q_a,'Time','Outflow','Q(t) vs t')

p1.grapher(t_a,T_res,'Time','Residence Time','t_res(t) vs t')

Part_C1()

def Part_C2():
r_top=np.arange(3,5.5,.5)
s=6.0/r_top
S_a=[]
for element in s:
Area=math.pi*(6.0/element)*math.sqrt(6.0**2+((6.0**2)/(element**2)))
S_a.append(Area)
pylab.figure(figsize=(8,6), dpi=80)
pylab.subplot(1,1,1)
pylab.plot(s, S_a, "ro:")
pylab.title('Surface Area vs s')
pylab.xlabel('s')
pylab.ylabel('Surface Area')
pylab.grid(which='major',axis='both')

Part_C2()

```

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.