%pylab inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from IPython.html.widgets import interact, interactive, fixed
Populating the interactive namespace from numpy and matplotlib
Here we will put a description of the code later.
%%file Toomre.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from IPython.html.widgets import interact, interactive
from scipy.optimize import fsolve
pi = np.pi
def ring(particles,radius,G,M):
'''
Set initial positions and velocities for all stars.
'''
particle = []
velocity = []
theta_n = 0
arclen = (2*pi)/particles ## Arc length equally divided amongst the number of particles around circle
v = np.sqrt((G*M)/radius) ## Velocity for central force to stay in a circular orbit
while len(particle) < particles:
angle = theta_n*arclen
beta = angle + (pi/2.) ## Angle beta = angle of the position plus 90 degrees, for velocity vector
theta_n += 1
particle.append((radius*np.cos(angle), radius*np.sin(angle))) ## x = r*cos(theta) y = r*sin(theta)
velocity.append((v*np.cos(beta), v*np.sin(beta))) ## vx = v*cos(beta) vy = v*sin(beta)
return np.array(particle),np.array(velocity) ## Returns two arrays, particle position and velocity.
def init_rings_123(G,M):
'''
Arrange stars into rings located at a distance that is a
certain percentage of Rmin. Rmin is the minimum distance
between the disruptor galaxy and the disrupted galaxy.
This function is only used on direct, retrograde, and light
mass disruptor cases.
'''
ring1,velocity1 = ring(12,.2,G,M) ## The positions of the stars are dependent on details found in the paper by Toomre et al.
ring2,velocity2 = ring(18,.3,G,M)
ring3,velocity3 = ring(24,.4,G,M)
ring4,velocity4 = ring(30,.5,G,M)
ring5,velocity5 = ring(36,.6,G,M)
rings = np.array([ring1,ring2,ring3,ring4,ring5])
velocity = np.array([velocity1,velocity2,velocity3,velocity4,velocity5])
return rings,velocity ## Returns arrays of both the positions and velocity.
def init_rings_4(G,M):
'''
Arrange stars into rings located at a distance that is a
certain percentage of Rmin. Rmin is the minimum distance
between the disruptor galaxy and the disrupted galaxy. This
function is only used on the heavy mass disruptor case.
'''
ring1,velocity1 = ring(12,.12,G,M) ## The positions of the stars are dependent on details found in the paper by Toomre et al.
ring2,velocity2 = ring(18,.18,G,M)
ring3,velocity3 = ring(24,.24,G,M)
ring4,velocity4 = ring(30,.30,G,M)
ring5,velocity5 = ring(36,.36,G,M)
rings = np.array([ring1,ring2,ring3,ring4,ring5])
velocity = np.array([velocity1,velocity2,velocity3,velocity4,velocity5])
return rings,velocity ## Returns arrays of both the positions and velocity.
def unpack_rings_vel(rings,velocity):
'''
Make 4 arrays that hold the information for the x and y star positions, and the
x and y star velocities.
'''
rx_points = [] ## x-coordinates of all stars
ry_points = [] ## y-coordinates of all stars
vrx = [] ## initial x velocity of all stars
vry = [] ## initial y velocity of all stars
for ring in rings:
for point in ring:
rx_points.append(point[0])
ry_points.append(point[1])
for ring in velocity:
for point in ring:
vrx.append(point[0])
vry.append(point[1])
return np.array(rx_points), np.array(ry_points), np.array(vrx), np.array(vry) ## Returns numpy arrays of values
def derivgalaxy(y,t,M,S):
'''
Function extracts information from the 4 arrays in unpack_rings_vel
function and plugs the values into the given differential equations
that describe the motion of the stars and the motion of the disruptor
galaxy. The output is a single array.
'''
G = 4.302e-3 #pc(M_solar)^-1 (km/s)^2
vRx = y[0]
vRy = y[2]
Rx = y[1]
Ry = y[3]
vrx = y[4]
vry = y[6]
rx = y[5]
ry = y[7]
R = np.sqrt(Rx**2+Ry**2)
delta_x = (Rx-rx)
delta_y = (Ry-ry)
dvrx_dt = -G * ((M/np.sqrt(rx**2. + ry**2.)**3.)*rx - (S/np.sqrt(delta_x**2.+delta_y**2.)**3.)*delta_x #Given differential equation describing the motion of the stars.
+ (S/np.sqrt(Rx**2.+Ry**2.)**3.)*Rx)
dvry_dt = -G * ((M/np.sqrt(rx**2. + ry**2.)**3.)*ry - (S/np.sqrt(delta_x**2.+delta_y**2.)**3.)*delta_y
+ (S/np.sqrt(Rx**2.+Ry**2.)**3.)*Ry)
dvRx_dt = -G * ((M+S)/(np.sqrt(Rx**2+Ry**2))**3)*Rx #Given differential equation describing the motion of the Disruptor Galaxy.
dvRy_dt = -G * ((M+S)/(np.sqrt(Rx**2+Ry**2))**3)*Ry
return np.array([dvRx_dt, vRx, dvRy_dt, vRy, dvrx_dt, vrx, dvry_dt, vry])
def Make_Master_Array(Case = 1,Rx0 = -8, Ry0 = -9,Initial_velocity_X = 0.85,Initial_velocity_Y = 0.65,t = 11., M=330., S=330., dt = 0.01):
'''
The function takes the single array from derivgalaxy and plugs it into
the odeint solver. The function then filters out the information associated
with the positions of the disruptor galaxy and stars at all timesteps between
0 and 20 with 0.0075 intervals. The output is a 2 dimension matrix where the
columns represent positions of the disruptor galaxy and all of the stars, and the
rows represent a particular time.
'''
G = 4.302e-3 #pc(M_solar)^-1 (km/s)^2\
if Case ==1 or Case == 2 or Case == 3:
rings,velocity = init_rings_123(G,M) ## Chooses which function to run according to the example chosen
elif Case == 4:
rings,velocity = init_rings_4(G,M)
rx0,ry0,vrx_0,vry_0 = unpack_rings_vel(rings,velocity) ## Converts values determined above to 1-D arrays
vRx_0 = Initial_velocity_X ## Initial velocity of disruptor galaxy in x direction
vRy_0 = Initial_velocity_Y ## Initial velocity of disruptor galaxy in y direction
ts = np.arange(0.,t+0.1,0.0075)
MasterArray = []
for n in range(len(rx0)): ## Runs for all 120 particles in initial condition vectors.
output = odeint(derivgalaxy, np.array([vRx_0,Rx0,vRy_0,Ry0,vrx_0[n],rx0[n],vry_0[n],ry0[n]]),
ts, args=(M, S)) ## Solve the differential equation for each time index
##and output the position values of the stars and disruptor galaxy.
rx = output[:,5]
ry = output[:,7]
if n == 0:
Rx = output[:,1]
Ry = output[:,3]
MasterArray.append(Rx)
MasterArray.append(Ry)
MasterArray.append(rx)
MasterArray.append(ry)
else:
MasterArray.append(rx)
MasterArray.append(ry)
return MasterArray
def Make_Plots(results,t, dt):
'''
Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.
This is the Direct Passage situation.
'''
index = int(t/dt)
plt.figure(figsize = (7, 7))
plt.grid()
plt.xlim(-10,7)
plt.ylim(-16,15)
plt.plot(results[0][:index], results[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results[2*i][index], results[2*i + 1][index], 'ro', label = "Stars")
plt.show()
def Make_Plots_Green_Star(results, t, dt, GreenStar):
index = int(t/dt)
plt.figure(figsize = (7, 7))
plt.grid()
plt.xlim(-10,7)
plt.ylim(-16,15)
plt.plot(results[0][:index], results[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results[2*i][index], results[2*i + 1][index], 'ro', label = "Stars")
for i in range(GreenStar, GreenStar+1):
plt.plot(results[2*i][index], results[2*i + 1][index], 'go', label = "Highlighted Star")
plt.show()
def Generate_Data(dataset = 'all', save = True):
'''
Calculate data for all of the stars and disruptor galaxy at all timesteps.
'''
if dataset == 'all':
results_A = Make_Master_Array(Case = 1, Rx0 = -8, Ry0 = -9, Initial_velocity_X = 0.85,Initial_velocity_Y = 0.65,t = 20, M=330., S=330., dt = 0.0075)
#Direct Passage
results_B = Make_Master_Array(Case = 2, Rx0 = -8, Ry0 = 9,Initial_velocity_X = 0.85,Initial_velocity_Y = -0.65,t = 20, M=330., S=330., dt = 0.0075)
#Retrograde Passage
results_C = Make_Master_Array(Case = 3,Rx0 = -8, Ry0 = -9, Initial_velocity_X = 0.85,Initial_velocity_Y = 0.65,t = 20, M=330., S=82.5, dt = 0.0075)
#Light Mass Disruptor
results_D = Make_Master_Array(Case = 4, Rx0 = -8, Ry0 = -9, Initial_velocity_X = 0.85,Initial_velocity_Y = 0.65,t = 20, M=82.5, S=330., dt = 0.0075)
#Heavy Mass Disruptor
if save == True:
np.save('Toomre_A.npy', results_A)
np.save('Toomre_B.npy', results_B)
np.save('Toomre_C.npy', results_C)
np.save('Toomre_D.npy', results_D)
Overwriting Toomre.py
import Toomre as Toomre
Toomre.Generate_Data()
ls
Volume in drive C has no label. Volume Serial Number is 0839-5207 Directory of C:\Users\masha\summer2014\GalaxyProject 08/05/2014 01:15 PM <DIR> . 08/05/2014 01:15 PM <DIR> .. 07/14/2014 11:10 AM 7 .gitignore 06/26/2014 01:52 PM 108,063 1stAttempt.ipynb 07/24/2014 04:18 PM 23,781 1WackyStarPlot.ipynb 06/26/2014 02:06 PM 733,591 2ndAttempt.ipynb 07/17/2014 01:57 PM 712,906 36StarPlots.ipynb 06/26/2014 02:12 PM 450,200 3rdAttempt.ipynb 07/24/2014 02:39 PM 8,520 AllAccelerations.ipynb 07/30/2014 11:10 AM 9,153 DisruptorPassage.ipynb 07/21/2014 01:36 PM 32,333 EatenStars.ipynb 07/31/2014 01:42 PM 189 EllipticalOrbitInvestigation.ipynb 07/30/2014 02:49 PM 597,575 EnergyPlots.ipynb 07/28/2014 11:54 AM 131,204 HeavyMass0.01.ipynb 07/25/2014 03:17 PM 5,975,390 HeavyMassInvestigation.ipynb 07/24/2014 04:39 PM 32,875 InvestigateAccelerations.ipynb 07/14/2014 11:11 AM 795 miscellaneous.ipynb 08/03/2014 08:34 PM 62,179 RevisedRevisedSpeedUpCode.ipynb 07/30/2014 12:00 PM 60,970 RevisedSpeedUpCode.ipynb 07/22/2014 03:46 PM 5,874,376 RhoVsT.ipynb 07/16/2014 11:32 AM 75,597 SpeedUpCode.ipynb 08/05/2014 01:15 PM 8,840 Toomre.py 08/05/2014 01:15 PM 7,267 Toomre.pyc 08/05/2014 01:17 PM 5,190,496 Toomre_A.npy 08/05/2014 01:17 PM 5,190,496 Toomre_B.npy 08/05/2014 01:17 PM 5,190,496 Toomre_C.npy 08/05/2014 01:17 PM 5,190,496 Toomre_D.npy 08/05/2014 01:15 PM 70,911 ToomreCode.ipynb 08/05/2014 01:13 PM 57,054 ToomreExamples.ipynb 07/22/2014 03:21 PM 35,713 Turn1StarGreen.ipynb 07/21/2014 12:29 PM 50,854 Unit_Investigation.ipynb 08/02/2014 04:01 PM 189 Untitled0.ipynb 07/24/2014 03:42 PM 5,156 WeeklyPlan07-14.ipynb 31 File(s) 35,887,672 bytes 2 Dir(s) 110,955,401,216 bytes free
test_array = np.load('Toomre_A.npy')
print test_array.shape
(242L, 2681L)
def Make_Plots(results,t, dt):
'''
Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.
This is the Direct Passage situation.
'''
index = int(t/dt)
plt.plot(results[0][:index], results[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results[2*i][index], results[2*i + 1][index], 'ro', label = "Stars")
plt.show()
#a = interact(Make_Plots,t = (0.,20.1), dt = (0.0075))
A = interact(Toomre.Make_Plots(results = np.load('Toomre_A.npy'), t = 20.1, dt = 0.0075), results = fixed(results), t = (0.,20.1), dt = fixed(0.0075))
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-2-ec77c5ddff45> in <module>() ----> 1 A = interact(Toomre.Make_Plots(results = np.load('Toomre_A.npy'), t = 20.1, dt = 0.0075), results = fixed(results), t = (0.,20.1), dt = fixed(0.0075)) NameError: name 'interact' is not defined
def Make_Plots_B(t, dt):
'''
Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.
This is the Retrograde Passage situation.
'''
index = int(t/dt)
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.plot(results_B[0][:index], results_B[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results_B[2*i][index], results_B[2*i + 1][index], 'ro', label = "Stars")
plt.show()
b = interact(Make_Plots_B,t = (0.,20.1), dt = (0.0075))
def Make_Plots_C(t, dt):
'''
Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.
This is the Light Mass Disruptor situation.
'''
index = int(t/dt)
plt.plot(results_C[0][:index], results_C[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results_C[2*i][index], results_C[2*i + 1][index], 'ro', label = "Stars")
plt.show()
c = interact(Make_Plots_C,t = (0.,20.1), dt = (0.0075))
def Make_Plots_D(t, dt):
'''
Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.
This is the Heavy Mass Disruptor situation.
'''
index = int(t/dt)
plt.plot(results_D[0][:index], results_D[1][:index], 'b--', label = 'Disturbant Galaxy')
for i in range(1,121):
plt.plot(results_D[2*i][index], results_D[2*i + 1][index], 'ro', label = "Stars")
plt.show()
d = interact(Make_Plots_D,t = (0.,20.1), dt = (0.0075))