import numpy as np import matplotlib.pyplot as plt import sys import random as random import mpl_toolkits.mplot3d.axes3d as p3 import matplotlib.animation as animation class part1: def __init__(self): self.H2gPrMol = 2.01588 # Hydrogen gass weight [g/mol] self.N_a = 6.02214179e23 # Avogadros number self.mass_sat = 1100 # Satellite/Rocket mass [kg] self.k = 1.380648e-23 #Boltzmann's constant [J/K] self.mass1Particle = self.H2gPrMol/self.N_a*1e-3 # [kg] def simulate_box(self, n, Nt, dt, box_dim, T): """ Calculates the motion of the particles used for fuel calculations variables: n = number of particles Nt = Number of time steps dt = Time step box_dim = Dimension of the box """ sigma = #... # Standard deviation for initial gaussian velocity distribution mean = 0 # Mean for initial gaussian distribution (no bulk flow) pos = np.zeros((n, 3, Nt)) # Position vector vel = np.zeros((n, 3)) # Velocity vector """ Setting initial positions and velocities. Positions are selected from a uniform distribution while velocities are selected from a gaussian distribution with a mean and deviation described in part1. """ for j in range(n): for k in range(3): pos[j, k, 0] = #... uniform dist. vel[j, k] = #... Gauss dist. x_momentum = 0.0 # Momentum value scalar in x-direction no_escapes = 0 # Counter for number of particles escaped for i in range(Nt-1): #Time loop pos[:,:,i+1] = pos[:,:,i] + vel[:,:]*dt #Move ALL particles one step #If you are really skilled with vectorization, #It's possible to vectorize the particle loop too (hint: masking) for j in range(n): # Particle loop """ Here you correct for collisions with the walls. You might need som if-else tests to make sure you cover all possible collisions """ #Potential solution """ #Check if particle is outside of any boundary if pos[j, 0, i+1] >= box_dim: # we are outside the box in x-direction if vel[j, 0] > 0: # the velocity is pointing outwards #Here turn the direction of the velocity (bounce particle) #This wall is special, it has a hole! Check if we hit it. if (((pos[j,1,i+1] > hole_min) and (pos[j,1,i+1] < hole_max)) # within hole in y-direction and (...) ): #particle is within hole in z-directin no_escapes += 1 #One particle escaped x_momentum += ... #Here calculate additional momentum else: continue #particle already bounced, don't turn it. elif pos[j,0, i+1] <= 0: #Similar test for negative x-dir if vel[j, 0] < 0: # the velocity is pointing outwards #Turn the direction of the velocity (bounce particle) else: continue #particle already bounced! elif pos[j,1, i+1] >= box_dim: # positive y-dir ... ... """ return x_momentum, pos def launch_sim(self, no_boxes, v_end, A, no_particles_sec): """ Simulates the launch. Variables: no_boxes: the number of fuel boxes you need v_end: the desired velocity (relative to surface) at the end of launch A: dp/dt for the box described in first method no_particles_sec: number of particles pr second leaving the box in first method returns: How much time it takes to achieve the boost and how much fuel you have used up. """ A = A*no_boxes # (dp/dt) for whole engine B = 0. #(insert value) # dm/dt for the box described in first method v = 0. # initial velocity relative to surface of planet time = 0. # initialize time variable T = 20.*60 # Total time, 20 minutes Nt = 10000 # Number of time steps dt = float(T)/Nt # Time step initial_fuel_mass = 100 # Calculate/set fuel mass M = self.mass_sat + initial_fuel_mass # Total mass for i in xrange(Nt): """ Here you need to update the new value of the velocity aswell as the new total mass value M """ #v += boost*dt - gravity ... #M -= ... time += dt if M < self.mass_sat: """ Some sort of error message telling you you're out of fuel """ elif v < 0: """ Gravity stronger than engine """ elif v >= v_end: """ Boost is successful. Save values and end method """ #fuel_needed = ... print "Boost Succesful" return fuel_needed, time return 0., 0. #returns 0 because the boost was not successful. def plot(self, n, Nt, dt, box_dim, T): """ A method for plotting and animating a simulation of particles in box """ def update_lines(num, dataLines, lines) : for line, data in zip(lines, dataLines) : line.set_data(data[0:2, num-1:num]) line.set_3d_properties(data[2,num-1:num]) return lines # Attach 3D axis to the figure fig = plt.figure() ax = p3.Axes3D(fig) m = 100 #Run the actual simulation x_momentum, datax = self.simulate_box(n, Nt, dt, box_dim, T) lines = [] data = [] for i in range(n): data.append( [datax[i]]) #wrap data inside another layer of [], needed for animation! lines.append([ax.plot(data[i][0][0,0:1], data[i][0][1,0:1], data[i][0][2,0:1], 'o')[0]]) # Set the axes properties ax.set_xlim3d([0.0, box_dim]) ax.set_xlabel('X') ax.set_ylim3d([0.0, box_dim]) ax.set_ylabel('Y') ax.set_zlim3d([0.0, box_dim]) ax.set_zlabel('Z') ax.set_title('Particle Animation') # Creating the Animation object ani = [i for i in range(n)] #"Initialize" a list of length n for i in range(n): #This is the method that needs the elements of data and lines to be wrapped in two layers of [] ani[i] = animation.FuncAnimation(fig, update_lines, m, fargs=(data[i], lines[i]), interval=50, blit=False) plt.show() instance = part1() instance.plot(box_dim = 1e-6, n = 100, Nt = 1000, dt = 1e-12, T = 10000)