{ "metadata": { "name": "", "signature": "sha256:a06f9ea2114e5d70ee39839c1ffad6f450c66ee6ec2fd9c9cfe980fd5a3c4ae8" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Toomre Code" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%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" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we will put a description of the code later." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%%file Toomre.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\n", "from scipy.optimize import fsolve\n", "\n", "pi = np.pi\n", "def ring(particles,radius,G,M):\n", " \n", " '''\n", " Set initial positions and velocities for all stars.\n", " '''\n", " \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))) ## vx = v*cos(beta) vy = v*sin(beta)\n", " return np.array(particle),np.array(velocity) ## Returns two arrays, particle position and velocity.\n", "\n", "def init_rings_123(G,M):\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.\n", " This function is only used on direct, retrograde, and light\n", " mass disruptor cases.\n", " '''\n", " \n", " ring1,velocity1 = ring(12,.2,G,M) ## The positions of the stars 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 the positions and velocity.\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", " \n", " '''\n", " Make 4 arrays that hold the information for the x and y star positions, and the\n", " x and y star velocities. \n", " '''\n", " \n", " rx_points = [] ## x-coordinates of all stars\n", " ry_points = [] ## y-coordinates of all stars\n", " vrx = [] ## initial x velocity of all stars\n", " vry = [] ## initial y velocity of all stars\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", " \n", " '''\n", " Function extracts information from the 4 arrays in unpack_rings_vel\n", " function and plugs the values into the given differential equations\n", " that describe the motion of the stars and the motion of the disruptor\n", " galaxy. The output is a single array.\n", " '''\n", " \n", " G = 4.302e-3 #pc(M_solar)^-1 (km/s)^2\n", " vRx = y[0]\n", " vRy = y[2]\n", " Rx = y[1]\n", " Ry = y[3]\n", " vrx = y[4]\n", " vry = y[6]\n", " rx = y[5]\n", " ry = y[7]\n", " R = np.sqrt(Rx**2+Ry**2)\n", " delta_x = (Rx-rx)\n", " delta_y = (Ry-ry)\n", " \n", " 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.\n", " + (S/np.sqrt(Rx**2.+Ry**2.)**3.)*Rx)\n", " dvry_dt = -G * ((M/np.sqrt(rx**2. + ry**2.)**3.)*ry - (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 #Given differential equation describing the motion of the Disruptor Galaxy.\n", " dvRy_dt = -G * ((M+S)/(np.sqrt(Rx**2+Ry**2))**3)*Ry\n", " \n", " return np.array([dvRx_dt, vRx, dvRy_dt, vRy, dvrx_dt, vrx, dvry_dt, vry])\n", "\n", "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):\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", " \n", " G = 4.302e-3 #pc(M_solar)^-1 (km/s)^2\\\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", " \n", " rx0,ry0,vrx_0,vry_0 = unpack_rings_vel(rings,velocity) ## Converts values determined above to 1-D arrays\n", " vRx_0 = Initial_velocity_X ## Initial velocity of disruptor galaxy in x direction\n", " vRy_0 = Initial_velocity_Y ## Initial velocity of disruptor galaxy in y direction\n", " \n", " ts = np.arange(0.,t+0.1,0.0075)\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([vRx_0,Rx0,vRy_0,Ry0,vrx_0[n],rx0[n],vry_0[n],ry0[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[:,5] \n", " ry = output[:,7]\n", " \n", " if n == 0:\n", " \n", " Rx = output[:,1] \n", " Ry = output[:,3] \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_Plots(results,t, dt):\n", " \n", " '''\n", " Function extracts all positions of stars and the disruptor galaxy from a matrix and plots them.\n", " This is the Direct Passage situation.\n", " '''\n", " \n", " index = int(t/dt)\n", " plt.figure(figsize = (7, 7))\n", " plt.grid()\n", " plt.xlim(-10,7)\n", " plt.ylim(-16,15)\n", " plt.plot(results[0][:index], results[1][:index], 'b--', label = 'Disturbant Galaxy')\n", " for i in range(1,121):\n", " plt.plot(results[2*i][index], results[2*i + 1][index], 'ro', label = \"Stars\")\n", " plt.show()\n", " \n", "def Make_Plots_Green_Star(results, t, dt, GreenStar):\n", " index = int(t/dt)\n", " plt.figure(figsize = (7, 7))\n", " plt.grid()\n", " plt.xlim(-10,7)\n", " plt.ylim(-16,15)\n", " plt.plot(results[0][:index], results[1][:index], 'b--', label = 'Disturbant Galaxy')\n", " for i in range(1,121):\n", " plt.plot(results[2*i][index], results[2*i + 1][index], 'ro', label = \"Stars\")\n", " for i in range(GreenStar, GreenStar+1):\n", " plt.plot(results[2*i][index], results[2*i + 1][index], 'go', label = \"Highlighted Star\")\n", " plt.show()\n", " \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 = -8, Ry0 = -9, Initial_velocity_X = 0.85,Initial_velocity_Y = 0.65,t = 20, M=330., S=330., dt = 0.0075)\n", " #Direct Passage \n", "\n", " 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)\n", " #Retrograde Passage\n", "\n", " 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)\n", " #Light Mass Disruptor\n", "\n", " 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)\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)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Overwriting Toomre.py\n" ] } ], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "import Toomre as Toomre\n", " " ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "Toomre.Generate_Data()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "ls" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " Volume in drive C has no label.\n", " Volume Serial Number is 0839-5207\n", "\n", " Directory of C:\\Users\\masha\\summer2014\\GalaxyProject\n", "\n", "08/05/2014 01:15 PM