from yt.mods import *
import numpy as np
import pdb
import matplotlib.pyplot as plt
import yt.units as units
from yt.units import kpc
from yt.frontends.sph.data_structures import ParticleDataset
ParticleDataset.filter_bbox = True
ParticleDataset._skip_cache = True
Load the snapshot on a coarse box and take a projection plot
fname = '/Volumes/pegasus/gadgetruns/m13_mr_Dec16_2013/snapdir_257/snapshot_257.0.hdf5'

unit_base = {'UnitLength_in_cm'         : 3.08568e+21,
             'UnitMass_in_g'            :   1.989e+43,
             'UnitVelocity_in_cm_per_s' :      100000}

bbox_lim = 1.e5

bbox = [[-bbox_lim,bbox_lim],
        [-bbox_lim,bbox_lim],
        [-bbox_lim,bbox_lim]]

ds1 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds1.index
ad1= ds1.all_data()


print 'left edge: ',ds1.domain_left_edge
print 'right edge: ',ds1.domain_right_edge
print 'center: ',ds1.domain_center
print 'total SFR is [in entire volume]: ',np.sum(ad1[("PartType0","StarFormationRate")])

px = ProjectionPlot(ds1, 'y', ('deposit', 'all_density'))
px.show()
Find the density peak to serve as a center for the zoomed box
density = ad1[("PartType0","density")]
wdens = np.where(density == np.max(density))
coordinates = ad1[("PartType0","Coordinates")]
maxdens_coordinates = coordinates[wdens]
center = maxdens_coordinates[0]


bbox_lim = ds1.quan(5.e3,'code_length')

bbox = [[center[0]-bbox_lim,center[0]+bbox_lim],
        [center[1]-bbox_lim,center[1]+bbox_lim],
        [center[2]-bbox_lim,center[2]+bbox_lim]]

print bbox
Load the zoomed snap, and ensure that the bulk of the SFR is still encapsulated
ds2 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds2.index
ds2.periodicity=(False,False,False)
ad2= ds2.all_data()
np.sum(ad2["PartType0","StarFormationRate"])

print ds2.domain_center
print ds2.domain_width
print 'right edge: ',ds2.domain_right_edge
print 'left edge: ',ds2.domain_left_edge

px = ProjectionPlot(ds2, 'x', ('deposit', 'PartType0_smoothed_density'))
px.show()
SAME STUFF BUT DIFFERENT SNAPSHOT; compressed for readability
fname = '/Volumes/pegasus/gadgetruns/m13_mr_Dec16_2013/snapdir_239/snapshot_239.0.hdf5'
unit_base = {'UnitLength_in_cm'         : 3.08568e+21,
             'UnitMass_in_g'            :   1.989e+43,
             'UnitVelocity_in_cm_per_s' :      100000}

bbox_lim = 1.e5

bbox = [[-bbox_lim,bbox_lim],
        [-bbox_lim,bbox_lim],
        [-bbox_lim,bbox_lim]]
 

ds1 = load(fname,unit_base=unit_base,bounding_box=bbox)
ds1.index
ad1= ds1.all_data()


density = ad1[("PartType0","density")]
wdens = np.where(density == np.max(density))
coordinates = ad1[("PartType0","Coordinates")]
maxdens_coordinates = coordinates[wdens]

center = maxdens_coordinates[0]

bbox_lim = ds1.quan(5.e3,'code_length')

bbox2 = [[center[0]-bbox_lim,center[0]+bbox_lim],
        [center[1]-bbox_lim,center[1]+bbox_lim],
        [center[2]-bbox_lim,center[2]+bbox_lim]]

print bbox2

ds2 = load(fname,unit_base=unit_base,bounding_box=bbox2)
ds2.index
ds2.periodicity=(False,False,False)
ad2= ds2.all_data()
print "zoomed SFR: ",np.sum(ad2["PartType0","StarFormationRate"])

print ds2.domain_center
print ds2.domain_width
print 'right edge: ',ds2.domain_right_edge
print 'left edge: ',ds2.domain_left_edge

px = ProjectionPlot(ds2, 'x', ('deposit', 'PartType0_smoothed_density'))
px.show()