{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "from IPython.html.widgets import interact, interactive, fixed\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting ToomreGalaxy.py\n" ] } ], "source": [ "%%file ToomreGalaxy.py\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "from IPython.html.widgets import interact, interactive, fixed\n", "\n", "G = 4.49955370898e-08 #kpc^3 Msol^-1 (10^8 years)^-2\n", "RCM = 0.; xCM = 0.; yCM = 0. #let origin be at CM\n", "Rmin = 25 #kpc\n", "pi = np.pi\n", "\n", "def Set_Init_R_Cond(Ry0, M, S):\n", " #minimum separation distance\n", " Rmin = 25.0 #kpc\n", "\n", " #Velocity at distance of closest approach\n", " Vmin = np.sqrt(2.*G*(M+S)/Rmin) #parabolic orbit\n", "\n", " #Angular momentum per unit mass at distance of closest approach\n", " hmin = Rmin*Vmin #r x v - angular momentum per unit mass \n", "\n", " #From the orbital equations, relationship between geometric and physical parameters\n", " c = hmin**2/(G*(M+S))\n", " beta = c/2.\n", " alpha = -1./(2*c)\n", "\n", " #Set a range of y values and compute the corresponding x,R values for initial points on a parabolic orbit\n", " \n", " Rx0 = beta + alpha*Ry0**2\n", "\n", " R0 = np.sqrt(Rx0**2+Ry0**2)\n", "\n", " #unit tangent vector for the parabola\n", " T_x = 2*alpha*Ry0/np.sqrt((2*alpha*Ry0)**2+1)\n", " T_y = 1./np.sqrt((2*alpha*Ry0)**2+1)\n", "\n", " #Velocity vector\n", " V0 = np.sqrt(2.*G*(M+S)/R0)\n", " Vx0 = V0*T_x\n", " Vy0 = V0*T_y\n", " \n", " if (Ry0>0.):\n", " Vx0 *= -1.\n", " Vy0 *= -1.\n", " return Rx0, Ry0, Vx0, Vy0\n", "\n", "def ring(particles,percent,G,M):\n", " \n", " radius = Rmin * percent\n", " particle = []\n", " velocity = []\n", " theta_n = 0\n", " arclen = (2*pi)/particles ## Arc length equally divided amongst the number of particles around circle\n", " v = np.sqrt((G*M)/radius) ## Velocity for central force to stay in a circular orbit\n", " while len(particle) < particles:\n", " angle = theta_n*arclen\n", " beta = angle + (pi/2.) ## Angle beta = angle of the position plus 90 degrees, for velocity vector\n", " theta_n += 1\n", " particle.append((radius*np.cos(angle), radius*np.sin(angle))) ## x = r*cos(theta) y = r*sin(theta)\n", " velocity.append((v*np.cos(beta), v*np.sin(beta))) ## Same concept here as above.\n", " return np.array(particle),np.array(velocity) ## Returns two arrays, particle position and velocity.\n", "\n", "def init_rings_123(G,M):\n", " ring1,velocity1 = ring(12,.2,G,M) ## All of these are dependent on details found in the paper by Toomre et al.\n", " ring2,velocity2 = ring(18,.3,G,M)\n", " ring3,velocity3 = ring(24,.4,G,M)\n", " ring4,velocity4 = ring(30,.5,G,M)\n", " ring5,velocity5 = ring(36,.6,G,M)\n", " rings = np.array([ring1,ring2,ring3,ring4,ring5])\n", " velocity = np.array([velocity1,velocity2,velocity3,velocity4,velocity5])\n", " return rings,velocity ## Returns arrays of both\n", "\n", "def init_rings_4(G,M):\n", " \n", " '''\n", " Arrange stars into rings located at a distance that is a \n", " certain percentage of Rmin. Rmin is the minimum distance \n", " between the disruptor galaxy and the disrupted galaxy. This \n", " function is only used on the heavy mass disruptor case.\n", " '''\n", " \n", " ring1,velocity1 = ring(12,.12,G,M) ## The positions of the stars are dependent on details found in the paper by Toomre et al.\n", " ring2,velocity2 = ring(18,.18,G,M)\n", " ring3,velocity3 = ring(24,.24,G,M)\n", " ring4,velocity4 = ring(30,.30,G,M)\n", " ring5,velocity5 = ring(36,.36,G,M)\n", " rings = np.array([ring1,ring2,ring3,ring4,ring5])\n", " velocity = np.array([velocity1,velocity2,velocity3,velocity4,velocity5])\n", " return rings,velocity ## Returns arrays of both the positions and velocity.\n", "\n", "def unpack_rings_vel(rings,velocity):\n", " rx_points = [] ## x-coordinates of all massless particles\n", " ry_points = [] ## y-coordinates\n", " vrx = [] ## initial x velocity\n", " vry = [] ## initial y velocity\n", " for ring in rings:\n", " for point in ring:\n", " rx_points.append(point[0])\n", " ry_points.append(point[1])\n", " for ring in velocity:\n", " for point in ring:\n", " vrx.append(point[0])\n", " vry.append(point[1])\n", " return np.array(rx_points), np.array(ry_points), np.array(vrx), np.array(vry) ## Returns numpy arrays of values\n", "\n", "def derivgalaxy(y,t,M,S):\n", " G = 4.49955370898e-08 #kpc^3 M_sol^-1 unitTime^-2\n", " Rx = y[0]\n", " Vx = y[1]\n", "\n", " Ry = y[2]\n", " Vy = y[3]\n", " \n", " rxs = y[4]\n", " vrx = y[5]\n", " \n", " rys = y[6]\n", " vry = y[7]\n", " \n", " delta_x = (Rx-rxs)\n", " delta_y = (Ry-rys)\n", " \n", " R = np.sqrt(Rx**2+Ry**2)\n", " \n", " dvrx_dt = -G * ((M/np.sqrt(rxs**2. + rys**2.)**3.)*rxs - (S/np.sqrt(delta_x**2.+delta_y**2.)**3.)*delta_x \n", " + (S/np.sqrt(Rx**2.+Ry**2.)**3.)*Rx)\n", " dvry_dt = -G * ((M/np.sqrt(rxs**2. + rys**2.)**3.)*rys - (S/np.sqrt(delta_x**2.+delta_y**2.)**3.)*delta_y \n", " + (S/np.sqrt(Rx**2.+Ry**2.)**3.)*Ry)\n", " \n", " dvRx_dt = -G * ((M+S)/(np.sqrt(Rx**2+Ry**2))**3)*Rx \n", " dvRy_dt = -G * ((M+S)/(np.sqrt(Rx**2+Ry**2))**3)*Ry\n", " \n", " return np.array([Vx, dvRx_dt, Vy, dvRy_dt, vrx, dvrx_dt, vry, dvry_dt])\n", "\n", "def Make_Master_Array(Case = 1, Rx0 = -39., Ry0 = -80.,M=1.0e11, S=1.0e11):\n", " \n", " '''\n", " The function takes the single array from derivgalaxy and plugs it into\n", " the odeint solver. The function then filters out the information associated\n", " with the positions of the disruptor galaxy and stars at all timesteps between\n", " 0 and 20 with 0.0075 intervals. The output is a 2 dimension matrix where the \n", " columns represent positions of the disruptor galaxy and all of the stars, and the\n", " rows represent a particular time.\n", " '''\n", " Rmin = 25 #kpc\n", " tmax = 20.\n", " dt = 0.007\n", " ts = np.arange(0.,tmax+dt/10.,dt)\n", " \n", " \n", " if Case ==1 or Case == 2 or Case == 3:\n", " rings,velocity = init_rings_123(G,M) ## Chooses which function to run according to the example chosen\n", " elif Case == 4:\n", " rings,velocity = init_rings_4(G,M)\n", " \n", " Rx0, Ry0, Vx0, Vy0 = Set_Init_R_Cond(Ry0=Ry0, M = M, S = S)\n", " rx0,ry0,vx0,vy0 = unpack_rings_vel(rings,velocity) ## Converts values determined above to 1-D arrays\n", " \n", " \n", " MasterArray = []\n", " \n", " for n in range(len(rx0)): ## Runs for all 120 particles in initial condition vectors.\n", " \n", " output = odeint(derivgalaxy, np.array([Rx0, Vx0, Ry0, Vy0, rx0[n], vx0[n],ry0[n], vy0[n]]),\n", " ts, args=(M, S)) ## Solve the differential equation for each time index \n", " ##and output the position values of the stars and disruptor galaxy.\n", " \n", " \n", " rx = output[:,4] \n", " ry = output[:,6]\n", " \n", " if n == 0:\n", " \n", " Rx = output[:,0] \n", " Ry = output[:,2] \n", " \n", " MasterArray.append(Rx)\n", " MasterArray.append(Ry)\n", " MasterArray.append(rx)\n", " MasterArray.append(ry)\n", " \n", " \n", " else:\n", " MasterArray.append(rx)\n", " MasterArray.append(ry)\n", " \n", " return MasterArray\n", "\n", "def Make_Plot_stars(results, M, S, t, dt):\n", "\n", " index = int(t/dt)\n", " \n", " Rx = results[0][:index]\n", " Ry = results[1][:index]\n", " RxS = xCM + (M/(M+S))*Rx\n", " RyS = yCM + (M/(M+S))*Ry\n", " RxM = xCM - (S/(M+S))*Rx\n", " RyM = yCM - (S/(M+S))*Ry\n", " RxS -= RxM\n", " RyS -= RyM\n", " RxM -= RxM\n", " RyM -= RyM\n", " plt.plot(RxS, RyS, 'b--', label = 'Disturbing Galaxy')\n", " plt.plot(RxS[-1], RyS[-1], 'bo')\n", " plt.plot(RxM, RyM, 'r--', label = 'Main Galaxy')\n", " plt.plot(RxM[-1], RyM[-1], 'ro')\n", " for i in range(1, 121):\n", " plt.plot(results[2*i][index]+RxM[-1], results[2*i + 1][index]+RyM[-1], 'go', label = \"Stars\")\n", " \n", " plt.xlim(1.1*Rx[0],xCM-1.1*Rx[0])\n", " plt.ylim(1.1*Ry[0],yCM-1.1*Ry[0])\n", " plt.xlim(-150,150)\n", " plt.ylim(-150,150)\n", " plt.grid()\n", " #plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),\n", " #ncol=2, fancybox=True, shadow=True)\n", " plt.show()\n", " \n", "def Make_Plots_Yellow_Star(results, M, S, t, dt, YellowStar):\n", " \n", " index = int(t/dt)\n", " \n", " Rx = results[0][:index]\n", " Ry = results[1][:index]\n", " RxS = xCM + (M/(M+S))*Rx\n", " RyS = yCM + (M/(M+S))*Ry\n", " RxM = xCM - (S/(M+S))*Rx\n", " RyM = yCM - (S/(M+S))*Ry\n", " RxS -= RxM\n", " RyS -= RyM\n", " RxM -= RxM\n", " RyM -= RyM\n", " plt.plot(RxS, RyS, 'b--', label = 'Disturbing Galaxy')\n", " plt.plot(RxS[-1], RyS[-1], 'bo')\n", " plt.plot(RxM, RyM, 'r--', label = 'Main Galaxy')\n", " plt.plot(RxM[-1], RyM[-1], 'ro')\n", " for i in range(1, 121):\n", " plt.plot(results[2*i][index]+RxM[-1], results[2*i + 1][index]+RyM[-1], 'go', label = \"Stars\")\n", " for i in range(YellowStar, YellowStar+1):\n", " plt.plot(results[2*i][index], results[2*i + 1][index], 'yo', label = \"Highlighted Star\") \n", " #plt.xlim(1.1*x[0],xCM-1.1*x[0])\n", " #plt.ylim(1.1*y[0],yCM-1.1*y[0])\n", " plt.xlim(-150,150)\n", " plt.ylim(-150,150)\n", " plt.grid()\n", " #plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),\n", " #ncol=2, fancybox=True, shadow=True)\n", " plt.show()\n", " \n", "def Generate_Data(dataset = 'all', save = True):\n", " \n", " '''\n", " Calculate data for all of the stars and disruptor galaxy at all timesteps.\n", " '''\n", " \n", " if dataset == 'all':\n", "\n", " results_A = Make_Master_Array(Case = 1, Rx0 = -39., Ry0 = -80., M=1.0e11, S=1.0e11) \n", " #Direct Passage \n", "\n", " results_B = Make_Master_Array(Case = 2, Rx0 = -39., Ry0 = 80., M=1.0e11, S=1.0e11) \n", " #Retrograde Passage\n", "\n", " results_C = Make_Master_Array(Case = 3, Rx0 = -39., Ry0 = -80., M=1.0e11, S=1.0e11/4)\n", " #Light Mass Disruptor\n", "\n", " results_D = Make_Master_Array(Case = 4, Rx0 = -39., Ry0 = -80., M=1.0e11, S=1.0e11*4)\n", " #Heavy Mass Disruptor\n", "\n", " if save == True:\n", " \n", " np.save('Toomre_A.npy', results_A)\n", " np.save('Toomre_B.npy', results_B)\n", " np.save('Toomre_C.npy', results_C)\n", " np.save('Toomre_D.npy', results_D) \n", " " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import ToomreGalaxy as ToomreGalaxy" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ToomreGalaxy.Generate_Data()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }