python/quick_check.py at main · hannahjleary/python · GitHub

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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import h5py
from mpl_toolkits.axes_grid1 import make_axes_locatable

mp = 1.672622e-24 # mass of hydrogren atom, in grams
kb = 1.380658e-16 # boltzmann constant in ergs/K
mu = 0.6 # mean molecular weight (mu) of 1

DE = 1 # Dual Energy Flag

dnamein='../../../../../ix/eschneider/hjl28/data/radiative/super/16retry2/hdf5/' # directory where the file is located
dnameout='../../../../../ix/eschneider/hjl28/plots/radiative/super/16retry2/'

CAT = 1

# t_cc = 4.89e4 # (vwind = 10 km/s)
# t_cc = 4.89e3 # cloud crushing time in kyr (vwind = 100 km/s)
t_cc = 4.89e2 # cloud crushing time in kyr (vwind = 1000 km/s)
istart = 0
iend = 500
step = 10
time = 0

for i in range(istart, iend, step):

    if CAT:
        f = h5py.File(dnamein + str(i) + '_slice.h5', 'r') # open the hdf5 file for reading
    else:
        f = h5py.File(dnamein + str(i) + '_slice.h5.0', 'r') # open the hdf5 file for reading
    head = f.attrs # read the header attributes into a structure, called head

    gamma = head['gamma'] # ratio of specific heats
    t  = head['t'] # time of this snapshot, in kyr
    nx = head['dims'][0] # number of cells in the x direction
    ny = head['dims'][1] # number of cells in the y direction
    nz = head['dims'][2] # number of cells in the z direction
    dx = head['dx'][0] # width of cell in x direction
    l_c = head['length_unit']
    t_c = head['time_unit']
    m_c = head['mass_unit']
    d_c = head['density_unit']
    v_c = head['velocity_unit']
    e_c = head['energy_unit']
    p_c = e_c # pressure units are the same as energy density units, density*velocity^2/length^3

    d  = f['d_xy'][:]
    px  = f['mx_xy'][:]
    py  = f['my_xy'][:]
    pz  = f['mz_xy'][:]
    E = f['E_xy'][:]

    if DE:
        GE = f['GE_xy'][:]

    f.close()

    n = d * d_c/ (mu*mp) # number density, particles per cm^3

    vx = px/d
    vy = py/d
    vz = pz/d

    if not DE:
        KE = 0.5 * d * (vx*vx + vy*vy + vz*vz)
        GE = E - KE

    T = GE*(gamma-1.0)*p_c / (n*kb) #temperature
    logT = np.log10(T)

    km = 1e-5

    Vx = vx*v_c*km #velocity in the x direction

    # if CAT:
    #     f = h5py.File(dnamein + str(i) + '_proj.h5', 'r') # open the hdf5 file for reading
    # else:
    #     f = h5py.File(dnamein + str(i) + '_proj.h5.0', 'r') # open the hdf5 file for reading
    # head = f.attrs # read the header attributes into a structure, called head
    # d  = f['d_xy'][:]
    # # T = f['T_xy'][:]


    # d = d * m_c / (l_c**2)
    # n = d / (mu*mp) # number density, particles per cm^3
    # logn = np.log10(n)


    # P = np.log10(n*kb*T)

    # f.close()

    # print(d_c)
    # print(p_c)
    # print('\t min \t\t\t max')
    # print('n: ', np.min(logn) , '\t' , np.max(logn))
    # print('T: ', np.min(logT) , '\t' , np.max(logT))
    # print('Vx: ', np.min(Vx) , '\t' , np.max(Vx))

    # mins = [Tmin, nmin, vmin]
    # maxs = [Tmax, nmax, vmax]
    cmaps = ['plasma', 'viridis', 'YlOrRd']
    labels = ['$log_{10}(K)$', '$log_{10}(N_{H})$ [$cm^{-2}$]', '$kms^{-1}$']

    fig = plt.figure()
    fig_color = 'white'
    bg_color = 'black'

    im = plt.imshow(logT.T, cmap=cmaps[0])
    plt.colorbar(im)

    plt.savefig(dnameout + str(i) + '.png', dpi=300,
                bbox_inches='tight', pad_inches = 0.2, facecolor=bg_color) #facecolor=bg_color
    plt.close(fig)