In [13]:
# Import Required Python Base Packages
import sys
import unittest
import numpy as np
import random
import collections
import os
import time as tp
import glob
# Create a Safety Mechanism for Not Having MatPlotLib Installed
try:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import Axes3D
HAS_MATPLOTLIB = True
matplotlib.rcParams['figure.figsize'] = (16, 6)
matplotlib.rcParams['font.size'] = (14)
except ImportError:
HAS_MATPLOTLIB = False
# Import the Amuse Base Packages
from amuse import datamodel
from amuse.units import nbody_system
from amuse.units import units
from amuse.units import constants
from amuse.datamodel import particle_attributes
from amuse.io import *
from amuse.lab import *
# Import the Amuse Stellar Packages
from amuse.ic.kingmodel import new_king_model
from amuse.ic.salpeter import new_salpeter_mass_distribution_nbody
# Import the Amuse Gravity & Close-Encounter Packages
from amuse.community.ph4.interface import ph4 as grav
from amuse.community.smalln.interface import SmallN
from amuse.community.kepler.interface import Kepler
from amuse.couple import multiples
# Importing cPickle/Pickle
try:
import cPickle as pickle
except:
import pickle
from tycho import analysis
In [4]:
def get_timestep_from_filenames(file_path1, file_path2):
file_name1 = (file_path1.split("/")[-1])[:-5]
file_name2 = (file_path2.split("/")[-1])[:-5]
time1 = float((file_name1.split("_t")[-1]))
time2 = float((file_name2.split("_t")[-1]))
dt = np.abs(time2-time1)
return dt
#get_timestep_from_filenames("/home/jglaser/Tycho/Runs/Mordin/160426/Mordin_t400.000000.hdf5", "/home/jglaser/Tycho/Runs/Mordin/160426/Mordin_t400.05000.hdf5")
In [5]:
def plot_cluster_2D(particles, step_number, current_time, cluster_name, converter, run_dir, show_plot=False, save_plot=True):
# Skip Over Saved Plots
#file_dir = "%s/Plots/2D_XY/" %(run_dir)
file_dir = "/home/draco/jthornton/Tycho/Plots/"+cluster_name+"/"
if not os.path.exists(file_dir):
os.makedirs(file_dir)
file_name = "%s_%08d.png" %(cluster_name, step_number)
if not os.path.isfile(file_dir+file_name):
# Split Objects into Stars & Planets
number_of_planets = len([x for x in particles if x.type == "planet"])
number_of_stars = len([x for x in particles if x.type == "star"])
stars = particles.sorted_by_attribute('type')[number_of_planets:].copy()
planets = particles.sorted_by_attribute('type')[:number_of_planets].copy()
# Get Create Plotting Variables
s_x = (converter.to_si(stars.x)).value_in(units.parsec)
s_y = (converter.to_si(stars.y)).value_in(units.parsec)
p_x = (converter.to_si(planets.x)).value_in(units.parsec)
p_y = (converter.to_si(planets.y)).value_in(units.parsec)
s_sizes = (converter.to_si(stars.mass)).value_in(units.MSun) * 30.0+5
p_sizes = (converter.to_si(planets.mass)).value_in(units.MJupiter) * 30.0+30
plot_time = (converter.to_si(current_time)).value_in(units.Myr)
# Make the Plot!
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['Tahoma']
matplotlib.rcParams['figure.figsize'] = (16,9)
subplot(111)#, axisbg='#000003')
#gca().set_aspect('equal', adjustable='box')
plt.axis([-25,25,-25,25])
#plt.axis([-2,2,-1.1,1.1])
plt.xticks(np.arange(-25, 25, 10))
plt.grid()
plt.scatter(p_x,p_y, marker='o', s=p_sizes, color='r')
plt.scatter(s_x,s_y, marker='*', s=s_sizes, color='k')
# Making the Labels
plt.tick_params(axis='both', which='major', labelsize=16)
plt.title('The %s Cluster (%i Stars, %i Planets)' %(cluster_name, number_of_stars, number_of_planets), fontsize=25)
plt.xlabel('X-Axis (pc)', fontsize=20)
plt.ylabel('Y-Axis (pc)', fontsize=20)
plt.text(-17.5, -20, '%.2f Myr' %(plot_time), style='italic',bbox={'facecolor':'blue', 'alpha':0.5, 'pad':15}, fontsize=15)
# Showing or Saving the Plot
if show_plot:
plt.show()
if save_plot:
plt.ioff()
plt.savefig(file_dir+file_name, format="png", dpi=150)
#plt.savefig(file_dir+file_name, format="png", dpi=75)
plt.clf()
plt.close('all')
In [49]:
def run_plotting(run_dir, cluster_name, converter, save_plot=True, show_plot=False):
"""
This is a function created to plot a 2D projection of the Cluster for every time-step.
The Function Returns: A bunch of plots that are either saved and/or displayed.
"""
print "[UPDATE] Starting Bulk Plotting (%s) ..." %(tp.strftime("%Y/%m/%d-%H:%M:%S", tp.gmtime()))
sys.stdout.flush()
# First get all of the times from the MasterParticleSet file names
search = glob.glob(run_dir+cluster_name+"*.hdf5")
i=0
current_time = []
for name in search:
time = name.split("_")
current_time.append(float(time[2][1:-5]) | nbody_system.time)
step_number = 0
for time in current_time:
file_name = "%s_MS_t%.3f.hdf5" %(cluster_name, time.number)
objects = read_set_from_file(run_dir+file_name, format="hdf5", close_file = True)
plot_cluster_2D(objects, step_number, time, cluster_name, converter, run_dir, show_plot, save_plot)
step_number+=1
print "[UPDATE] Finished Bulk Plotting (%s)!" %(tp.strftime("%Y/%m/%d-%H:%M:%S", tp.gmtime()))
sys.stdout.flush()
#run_plotting('/home/jglaser/Tycho/Runs/Garrus/160426/', 'Garrus', 0.05 | nbody.time)
In [7]:
def Import_Cluster(file_path, cluster_name):
"""
This is a function created to read in a Cluster's Stellar Parcticle Set.
The Function Returns: The Set of Stars, AMUSE Nbody Converter, Cluster Name, & King's Model W0 Value
"""
# Parses Variables from File_Name
file_name = (file_path.split("/")[-1])[:-5]
print file_name
file_format = ("hdf5")
w = 2.5
# N = int((file_name.split("_")[-2])[1:])
# total_mass = float((file_name.split("_")[-1])[1:]) | units.MSun
# Reads in the Set of Stars from the Specified File Path
stars = read_set_from_file(file_path, format=file_format, close_file=True)
number_of_stars = len(stars)
# This converter is set so that [1 mass = Cluster's Total Mass] & [1 length = 1 parsec].
# converter = nbody_system.nbody_to_si(total_mass, 1 | units.parsec)
return stars, cluster_name, w
In [8]:
def read_initial_state(cluster_name):
''' Reads in an initial state for the Tycho Module.
file_prefix: String Value for a Prefix to the Saved File
'''
# TODO: Convert the saved datasets from SI to NBody. Also everything else in this function.
# First, Define the Directory where Initial State is Stored
# file_dir = os.getcwd()+"/InitialState"
file_dir = "/home/draco/jthornton/Tycho/InitialState/"
file_base = file_dir + cluster_name
file_pkl = file_base + "_ic.pkl"
file_hdf5 = file_base + "_particles.hdf5"
# Second, Read the Master AMUSE Particle Set from a HDF5 File
file_format = "hdf5"
master_set = read_set_from_file(file_hdf5, format=file_format, close_file=True)
# Third, unPickle the Initial Conditions Array
ic_file = open(file_pkl, "rb")
ic_array = pickle.load(ic_file)
ic_file.close()
# Fourth convert ic_array.total_smass and viral_radius from strings to floats using string split
total_smass = float(ic_array.total_smass) | units.kg
viral_radius = float(ic_array.viral_radius) | units.m
# Fifth, Define the Master Set's Converter
converter = nbody_system.nbody_to_si(total_smass, viral_radius)
return master_set, ic_array, converter
In [38]:
w = 2.5
cluster_names = ["TestCluster4"]
runs_dir = "/home/draco/jthornton/Tycho/Run/MasterParticleSet/"
In [35]:
print runs_dir
print cluster_name
In [48]:
search = glob.glob(runs_dir+cluster_name+"*.hdf5")
i=0
current_time = []
for name in search:
time = name.split("_")
current_time.append(float(time[2][1:-5]))
print current_time
In [46]:
Energy, Time, T, U, L, P = analysis.GetValues("BinaryParticles2", start_time = | nbody_system.time, step_time = 0.125 | nbody_system.time)
analysis.EnergyGraph(Time, Energy, T, U, "BinaryStabilityTest2")
In [51]:
for cluster_name in cluster_names:
master_set, ic_array, converter = read_initial_state(cluster_name)
time_step = 0.03125
run_plotting(runs_dir, cluster_name, converter)
In [ ]:
Quick_Cluster = read_set_from_file("/home/draco/jthornton/BinaryTycho/Runs/400AU0Binaries2Run/160707/Master/400AU0Binaries2Run_Master_t0.050.hdf5", format="hdf5", close_file=True)
In [ ]:
plot_cluster_2D(Quick_Cluster, 5, 0.250 | nbody.time, "400AU0Binaries2Run", "/home/draco/jthornton/BinaryTycho/Runs/400AU0Binaries2Run/160707/", save_plot=False, show_plot=True)
In [60]:
particles = read_set_from_file("/home/draco/jthornton/Tycho/Run/MasterParticleSet/PlanetCluster2_MS_t0.125.hdf5", format="hdf5", close_file=True)
In [61]:
thing = particles.sorted_by_attribute('type')
print thing
In [ ]:
initial_state = "/home/draco/jthornton/Tycho/InitialState/movieTest2_particles.hdf5"
final_state = "/home/draco/jthornton/Tycho/Runs/MasterParticleSet/movieTest2_MS_t99.800.hdf5"
MasterSet1 = []
MasterSet2 = []
MasterSet1 = read_set_from_file(format = 'hdf5')
MasterSet2 = read_set_from_file()