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:
        #Ask Proff about this case
        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.

Delete this entry (admin only).