In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
gamma = 4.4983169634398596e-06
Defining a function that would determine my initial velocities of my star, $m$
The star's path is of circular nature, and can be explained with the equation for a circle centered at the origin: $x^{2}+y^2{2}=R^{2}$
Solving for $x$ gives: $x=\sqrt{R^{2}-y^{2}}$
The slope is then the derivative of this function:
$\frac{dx}{dy}=\text{-}y(R^{2}-y^{2})^{-\frac{1}{2}}$
This slope will serve to find the angle to split the magnitude of the velocity into $x$ and $y$ components
The star's initial velocity magnitude is determined by the equation for the velocity pertaining to a circular obrit: \begin{equation*} v = \sqrt{\frac{GM}{r}} \end{equation*} where $G$ is the gravitational constant, $M$ is the mass of the object orbited around, and $r$ is the distance from $m$ to $M$.
In [3]:
def velocities_m(M,x,y):
"""Computes the velocities in the x and y direction of a star, m
Parameters
-------------
M : central mass of main galaxy
x : position of star in x-direction
y : position of star in y-direction
Returns
------------
vx : velocity of star in x-direction
vy : velocity of star in y-direction
"""
rb = np.sqrt(x**2+y**2)
v = np.sqrt((gamma*M)/rb)
if x > 1e-5 or x < -1e-5:
dx_dy = -y*((rb**2-y**2)**-0.5)
else:
dx_dy = 0
if x < -1e-5 and y > 0 or x > 1e-5 and y < 0:
theta = np.arctan(abs(dx_dy))
elif x < -1e-5 and y < 0 or x > 1e-5 and y > 0:
theta = np.arctan(abs(1/dx_dy))
else:
theta = 0
if x < 0 and y > 0:
vx = -v*np.sin(theta)
vy = -v*np.cos(theta)
elif x < 0 and y < 0:
vx = v*np.cos(theta)
vy = -v*np.sin(theta)
elif x > 0 and y < 0:
vx = v*np.sin(theta)
vy = v*np.cos(theta)
elif x > 0 and y > 0:
vx = -v*np.cos(theta)
vy = v*np.sin(theta)
elif x == 0 and y > 0:
vx = -v
vy = 0
elif x == 0 and y < 0:
vx = v
vy = 0
elif x < 0 and y == 0:
vx = 0
vy = -v
elif x > 0 and y == 0:
vx = 0
vy = v
return vx,vy
Next, I define a function to compute the initial velocities of the disrupting galaxy, $S$
The disrupting galaxy S, follows a parabolic path around the main galaxy, $M$
The equation determining the shape of the parabola that $S$ will follow around M is: $x=-ay^{2}+by+c$
where $a=\frac{1}{4p}$, $b=\frac{-k}{2p}$, and $c=\frac{k^2}{4p}+h$
k,h, and p are as follows, where p is the distance from the vertex, (h,k), to the origin
In [5]:
y = np.linspace(-150,150,100)
plt.plot(-.01*y**2+25,y)
plt.scatter(0,0,color='y',label='origin')
plt.scatter(25,0,color='g',label='(h,k)')
plt.legend()
plt.ylim(-60,60)
plt.xlim(-25,35);
Assuming galaxy $M$'s frame of reference, $M$ is at the origin and $S$ follows the parabolic path around $M$
In this case then, $k=0$ and $h=p$, and $p$ is given in the original paper($25 \hspace{1 mm} kpc$), so the equation for the shape of the parabola that $S$ will follow simplifies to:
$x=\text{-}\frac{1}{4p}y^{2}+p$, which gives an easy way to figure out initial positions for galaxy S
The slope is then the derivative of this function:
$\frac{dx}{dy}=\text{-}\frac{y}{2p}$
This slope will serve to find the angle to split the magnitude of the velocity into $x$ and $y$ components
The magnitude of the velocity is determined by the equation pertaining to a parabolic orbit
$v = \sqrt{\frac{2\mu}{r}}$
where $\mu$ is the standard graitational parameter, $GM+GS$, and $r$ is the distance between the galaxies.
In [6]:
def velocities_S(M,S,x,y):
"""Computes the velocities in the x and y direction of a disrupting galaxy, S, following a parabolic orbit
Parameters
--------------
M : mass of main galaxy
S : mass of disrupting galaxy
x : position of S in x-direction
y : position of S in y-direction
Returns
--------------
vx : velocity of S in x-direction
vy : velocity of S in y-direction
"""
Rb = np.sqrt(x**2+y**2)
v = np.sqrt((2*gamma*(M+S))/(Rb))
dx_dy = (-y/50)
theta = np.arctan(abs(1/dx_dy))
if y > 0:
vx = v*np.cos(theta)
vy = -v*np.sin(theta)
elif y < 0:
vx = v*np.cos(theta)
vy = v*np.sin(theta)
return vx,vy