Paste #234

Welcome To LodgeIt

Welcome to the LodgeIt pastebin. In order to use the notification feature a 31 day cookie with an unique ID was created for you. The lodgeit database does not store any information about you, it's just used for an advanced pastebin experience :-). Read more on the about lodgeit page. Have fun :-)

hide this notification

Paste Details

reply | raw

posted on Jun 18, 2020, 8:13:34 PM

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import yt
from yt.utilities.cosmology import Cosmology
import matplotlib
matplotlib.use('Agg')
matplotlib.rc("font", size=18, family="serif")
from yt.mods import *
import sys
import numpy as np
from yt.units import mp
import yt.units as u
import matplotlib.pyplot as plt
from yt.utilities.physical_constants import *
from yt.data_objects.particle_filters import add_particle_filter

#Change absolute paths to files here
fbase = "group_subh_hgas_iso_hdf5_plt_cnt_0180"
fbase_part = "group_subh_hgas_iso_hdf5_part_0060"
out_dir = "./0060/"

def passive_part(pfilter, data):
    filter = data["all", "particle_mass"] < 0
    return filter

add_particle_filter("passive_part", function=passive_part, filtered_type='all',
                    requires=["particle_typepart"])

ds = load(fbase, particle_filename = fbase_part)

ds.add_particle_filter("passive_part")

def _passive_mass(field, data):
    return (data["particle_mass"]*-1.*1.E6*1.98E33)

ds.add_field("passive_mass", function=_passive_mass, units=r"g",
             display_name=r"Passive Mass", particle_type = True)    

#Field that stores cloud-in-cell mapped density of passive particles
def _passivecic_field(field, data):
    ptype = "passive_part"
    pos = data[ptype, 'particle_position']
    # Get back into density
    pden = data[ptype, 'passive_mass']#*1.E8*mass_sun_cgs
    d = data.deposit(pos, [pden], method = "cic")
    d = data.apply_units(d, "g")
    d /= data["index", "cell_volume"]
    return d

ds.add_field("passivecic_field", function=_passivecic_field, units=r"g*cm**-3",
             validators = [ValidateSpatial(0)], display_name=r"CIC Passive Field")

#projection
p = ProjectionPlot(ds, 'z', 'passivecic_field',center=[0., 0., 0.])
p.zoom(3)
p.set_log('passivecic_field', True)
p.save(out_dir+"0060_passivedens_z_proj.png")