diff --git a/astro_tigress/__init__.py b/astro_tigress/__init__.py index 5fc6f8a..2ebedd2 100644 --- a/astro_tigress/__init__.py +++ b/astro_tigress/__init__.py @@ -3,4 +3,6 @@ __license__ = "MIT" __description__ = "Python packages for the TIGRESS public data release" -from astro_tigress.tigress_read import * +from astro_tigress.tigress_read import Model, DataMHD, DataRad, DataChem, DataCO + +__all__ = ["Model", "DataMHD", "DataRad", "DataChem", "DataCO"] diff --git a/astro_tigress/const.py b/astro_tigress/const.py index 46781e0..85b81cf 100644 --- a/astro_tigress/const.py +++ b/astro_tigress/const.py @@ -11,225 +11,239 @@ # CGS PHYSICAL CONSTANTS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& -c = 2.99792458e10 # speed of light CGS -h = 6.6260755e-27 # Planck's constant CGS -g = 6.67259e-8 # Grav const CGS -kb = 1.380658e-16 # Boltzmann's const CGS -a = 7.56591e-15 # Radiation constant CGS -sb = 5.67051e-5 # sigma (stefan-boltzmann const) CGS -qe = 4.803206e-10 # Charge of electron CGS -ev = 1.60217733e-12 # Electron volt CGS -na = 6.0221367e23 # Avagadro's Number -me = 9.1093897e-28 # electron mass CGS -mp = 1.6726231e-24 # proton mass CGS -mn = 1.674929e-24 # neutron mass CGS -mh = 1.673534e-24 # hydrogen mass CGS -amu = 1.6605402e-24 # atomic mass unit CGS +c = 2.99792458e10 # speed of light CGS +h = 6.6260755e-27 # Planck's constant CGS +g = 6.67259e-8 # Grav const CGS +kb = 1.380658e-16 # Boltzmann's const CGS +a = 7.56591e-15 # Radiation constant CGS +sb = 5.67051e-5 # sigma (stefan-boltzmann const) CGS +qe = 4.803206e-10 # Charge of electron CGS +ev = 1.60217733e-12 # Electron volt CGS +na = 6.0221367e23 # Avagadro's Number +me = 9.1093897e-28 # electron mass CGS +mp = 1.6726231e-24 # proton mass CGS +mn = 1.674929e-24 # neutron mass CGS +mh = 1.673534e-24 # hydrogen mass CGS +amu = 1.6605402e-24 # atomic mass unit CGS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # ASTRONOMICAL CONSTANTS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # GENERAL -au = 1.496e13 # astronomical unit CGS -pc = 3.0857e18 # parsec CGS -yr = 3.155815e7 # sidereal year CGS -Myr = 1.0e6 * yr # million years CGS -ms = 1.98900e+33 # solar mass CGS -mj = 1.8986e30 # jupiter mass CGS -aj = 5.202545 * au # jupiter semi-major axis CGS -rs = 6.9599e10 # sun's radius CGS -ls = 3.839e33 # sun's luminosity CGS -mm = 7.35000e+25 # moon mass CGS -mer = 5.97400e+27 # earth mass CGS -rer = 6.378e8 # earth's radius CGS -mmsn_er = 1700.0 # g/cm^2, MMSN at 1AU -medd = 3.60271e+34 # Eddington mass CGS -km = 1.0e5 # km in cm +au = 1.496e13 # astronomical unit CGS +pc = 3.0857e18 # parsec CGS +yr = 3.155815e7 # sidereal year CGS +Myr = 1.0e6 * yr # million years CGS +ms = 1.98900e33 # solar mass CGS +mj = 1.8986e30 # jupiter mass CGS +aj = 5.202545 * au # jupiter semi-major axis CGS +rs = 6.9599e10 # sun's radius CGS +ls = 3.839e33 # sun's luminosity CGS +mm = 7.35000e25 # moon mass CGS +mer = 5.97400e27 # earth mass CGS +rer = 6.378e8 # earth's radius CGS +mmsn_er = 1700.0 # g/cm^2, MMSN at 1AU +medd = 3.60271e34 # Eddington mass CGS +km = 1.0e5 # km in cm # RADIO SPECIFIC -jy = 1.e-23 # Jansky CGS -restfreq_hi = 1420405751.786 # 21cm transition (Hz) -restfreq_co = 115271201800. # CO J=1-0 (Hz) -cm2perkkms_hi = 1.823e18 # HI column per intensity (thin) +jy = 1.0e-23 # Jansky CGS +restfreq_hi = 1420405751.786 # 21cm transition (Hz) +restfreq_co = 115271201800.0 # CO J=1-0 (Hz) +cm2perkkms_hi = 1.823e18 # HI column per intensity (thin) # OTHER WAVELENGTHS -ksun = 3.28 # abs K mag of sun (Binney & Merrifield) +ksun = 3.28 # abs K mag of sun (Binney & Merrifield) # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # GEOMETRY # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& -srdeg = 3283.0 # Degrees squared in a steradian -dtor = 0.0174532925199 # Degrees per radian -pi = 3.14159265359 # Pi +srdeg = 3283.0 # Degrees squared in a steradian +dtor = 0.0174532925199 # Degrees per radian +pi = 3.14159265359 # Pi # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # TELESCOPE STUFF # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& + def hpbw(lam_cm=None, freq_hz=None, diam_m=None, units=True): """ Return the half-power beam width for a telescope. """ - if (lam_cm == None and freq_hz == None) or \ - diam_m == None: + if (lam_cm is None and freq_hz is None) or diam_m is None: hpbw = None else: - if lam_cm == None: - lam_cm = c/freq_hz - hpbw = 1.2 * lam_cm / (diam_m*1e2) + if lam_cm is None: + lam_cm = c / freq_hz + hpbw = 1.2 * lam_cm / (diam_m * 1e2) - if units == True: - return {"units":"sr", - "value":hpbw} + if units: + return {"units": "sr", "value": hpbw} else: return hpbw + # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # PLANCK FUNCTION CALCULATIONS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& + def bb_inu(temp_k=None, lam_cm=None, units=True): """ Planck function. Returns f_nu given T and lambda. """ nu = c / lam_cm - fnu = 2.0*h*nu**3/c**2 /(np.exp(h*nu/(kb*temp_k))-1.0) + fnu = 2.0 * h * nu**3 / c**2 / (np.exp(h * nu / (kb * temp_k)) - 1.0) - if units == True: - return {"units":"erg/s/cm^2/sr/Hz", - "value":fnu} + if units: + return {"units": "erg/s/cm^2/sr/Hz", "value": fnu} else: return fnu + def bb_ilam(temp_k=None, lam_cm=None, units=True): """ Planck function. Returns f_lambda given T and lambda. """ - flam = 2.0*h*c**2/lam_cm**5 / (np.exp(h*c/(lam_cm*kb*temp_k))-1.0) + flam = 2.0 * h * c**2 / lam_cm**5 / (np.exp(h * c / (lam_cm * kb * temp_k)) - 1.0) - if units == True: - return {"units":"erg/s/cm^2/sr/cm", - "value":flam} + if units: + return {"units": "erg/s/cm^2/sr/cm", "value": flam} else: return flam + ### NOT FINISHED - FIX ### + def peak_inu(temp_k=None, units=True): """ Wavelength for peak of F_nu given T. NOT FUNCTIONAL YET. """ - peak_lam = 2.821439372/kb/temp_k/h - peak_nu = c / peak_lam + peak_lam = 2.821439372 / kb / temp_k / h + # peak_nu = c / peak_lam - if units == True: - return {"units":"cm", - "value":peak_lam} + if units: + return {"units": "cm", "value": peak_lam} else: return peak_lam + ### NOT FINISHED - FIX ### + def peak_ilam(temp_k=None, units=True): """ Wavelength for peak of F_lambda given T. NOT FUNCTIONAL YET. """ - peak_lam = h*c/kb/temp_k/4.965114231 - if units == True: - return {"units":"cm", - "value":peak_lam} + peak_lam = h * c / kb / temp_k / 4.965114231 + if units: + return {"units": "cm", "value": peak_lam} else: return peak_lam + # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # TRANSLATE HMS/DMS STRINGS BAND AND FORTH FROM DEC. DEG # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # CASA's QA tool should do this but sucks + def hms(val): """ Decimal value to 3-element tuple of hour, min, sec. Just calls dms/15 """ - return dms(val/15.0) + return dms(val / 15.0) + def dms(val): """ Decimal value to 3-element tuple of deg, min, sec. """ - if type(val) == type(np.array([1,2,3])): + if isinstance(val, type(np.array([1, 2, 3]))): val_abs = np.abs(val) d = np.floor(val_abs) - m = np.floor((val_abs - d)*60.0) - s = ((val_abs - d - m/60.0))*3600.0 - d *= +2.0*(val >= 0.0) - 1.0 + m = np.floor((val_abs - d) * 60.0) + s = (val_abs - d - m / 60.0) * 3600.0 + d *= +2.0 * (val >= 0.0) - 1.0 else: val_abs = abs(val) d = math.floor(val_abs) - m = math.floor((val_abs - d)*60.0) - s = ((val_abs - d - m/60.0))*3600.0 - d *= +2.0*(val >= 0.0) - 1.0 + m = math.floor((val_abs - d) * 60.0) + s = (val_abs - d - m / 60.0) * 3600.0 + d *= +2.0 * (val >= 0.0) - 1.0 return (d, m, s) + def ten(vec, hms=False): """ Convert 3-element vector to decimal degrees. Use HMS for hours. """ - if type(val) == type(np.array([1,2,3])): - dec = (np.abs(vec[0])+vec[1]/60.+vec[2]/3600.) - dec *= +2.0*(vec[0] >= 0.0) - 1.0 + if isinstance(vec, type(np.array([1, 2, 3]))): + dec = np.abs(vec[0]) + vec[1] / 60.0 + vec[2] / 3600.0 + dec *= +2.0 * (vec[0] >= 0.0) - 1.0 else: - dec = (abs(vec[0])+vec[1]/60.+vec[2]/3600.) - dec *= +2.0*(vec[0] >= 0.0) - 1.0 - if hms == True: - return 15.0*dec + dec = abs(vec[0]) + vec[1] / 60.0 + vec[2] / 3600.0 + dec *= +2.0 * (vec[0] >= 0.0) - 1.0 + if hms: + return 15.0 * dec else: return dec -def string_sixty(inp_str,hms=False): + +def string_sixty(inp_str, hms=False): """ Convert an HMS/DMS style string to a numeric vector. Again, QA should do this. """ - if hms==True: + if hms: pass else: pass -def sixty_string(inp_val,hms=False,colons=False): + +def sixty_string(inp_val, hms=False, colons=False): """ Convert a numeric vector to a string. Again, QA should do this. """ - if hms==True: + if hms: out_str = "%02dh%02dm%05.3fs" % (inp_val[0], inp_val[1], inp_val[2]) - elif colons==True: + elif colons: out_str = "%+02d:%02d:%05.3f" % (inp_val[0], inp_val[1], inp_val[2]) else: out_str = "%+02dd%02dm%05.3fs" % (inp_val[0], inp_val[1], inp_val[2]) return out_str + # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # LAPLACE COEFFICIENT AND ITS DIRIVATIVES CALCULATIONS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& + def binte(phi, n, s, alpha): """the integration of Laplace coefficient""" - return np.cos(n*phi) / (1 - 2*alpha*np.cos(phi) + alpha**2)**s + return np.cos(n * phi) / (1 - 2 * alpha * np.cos(phi) + alpha**2) ** s + def laplace_b(n, s, alpha): """the Laplace coefficient b^n_s(alpha)""" - return scipy.integrate.quad(binte, 0., 2.*math.pi, (n, s, alpha) - )[0] / math.pi + return scipy.integrate.quad(binte, 0.0, 2.0 * math.pi, (n, s, alpha))[0] / math.pi + + def dlaplace_b_dalpha(n, s, alpha): def deri(alpha0, n0, s0): return laplace_b(n0, s0, alpha0) - return scipy.misc.derivative(deri, alpha, dx=1.e-5, args=(n, s)) + + return scipy.misc.derivative(deri, alpha, dx=1.0e-5, args=(n, s)) + # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& # KEPLERIAN ORBIT COEFFICIENTS IN CGS # %&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%&%& def nkep(a, mstar=ms): """mean motion of orbit at semi-major axis a in CGS""" - return math.sqrt(g * mstar / float(a)**3) + return math.sqrt(g * mstar / float(a) ** 3) diff --git a/astro_tigress/radmc.py b/astro_tigress/radmc.py index b7caa38..335a37f 100644 --- a/astro_tigress/radmc.py +++ b/astro_tigress/radmc.py @@ -7,21 +7,20 @@ import yt import yt.units as u from astropy.io import fits -import copy from . import const from . import ytathena as ya from . import map_circular_beam as mcb + class YTData: def __init__(self, filename, left_edge=None, dims=None): self.fn = filename - fn_ext = filename.split('.')[-1] + fn_ext = filename.split(".")[-1] if fn_ext == "athdf": self.ytds = ya.load_athenapp(filename) else: - raise RuntimeError("file extension not recogonized: "+ - fn_ext) + raise RuntimeError("file extension not recogonized: " + fn_ext) if left_edge is None: self.left_edge = self.ytds.domain_left_edge else: @@ -30,19 +29,24 @@ def __init__(self, filename, left_edge=None, dims=None): self.dims = self.ytds.domain_dimensions else: self.dims = dims - self.ytdata = self.ytds.covering_grid(level=0, left_edge=self.left_edge, - dims=self.dims) + self.ytdata = self.ytds.covering_grid( + level=0, left_edge=self.left_edge, dims=self.dims + ) self.shape_xyz = self.dims self.nx = self.shape_xyz[0] self.ny = self.shape_xyz[1] self.nz = self.shape_xyz[2] self.dx = (self.ytdata["dx"][0, 0, 0].in_cgs()).value - nH0 = (self.ytdata["density"]/2.38858753789e-24/u.g*u.cm**3).in_cgs().to_ndarray() + nH0 = ( + (self.ytdata["density"] / 2.38858753789e-24 / u.g * u.cm**3) + .in_cgs() + .to_ndarray() + ) self.nH = np.swapaxes(nH0, 0, 2) - #velocity in km/s - velx0 = (self.ytdata["velocity_x"] / (u.km/u.s)).in_cgs().to_ndarray() - vely0 = (self.ytdata["velocity_y"] / (u.km/u.s)).in_cgs().to_ndarray() - velz0 = (self.ytdata["velocity_z"] / (u.km/u.s)).in_cgs().to_ndarray() + # velocity in km/s + velx0 = (self.ytdata["velocity_x"] / (u.km / u.s)).in_cgs().to_ndarray() + vely0 = (self.ytdata["velocity_y"] / (u.km / u.s)).in_cgs().to_ndarray() + velz0 = (self.ytdata["velocity_z"] / (u.km / u.s)).in_cgs().to_ndarray() vel = np.zeros([self.nx, self.ny, self.nz, 3]) for ix in range(self.nx): for iy in range(self.ny): @@ -51,31 +55,36 @@ def __init__(self, filename, left_edge=None, dims=None): vel[ix, iy, iz, 1] = vely0[ix, iy, iz] vel[ix, iy, iz, 2] = velz0[ix, iy, iz] self.vel_kms = np.swapaxes(vel, 0, 2) - #abd: CO, H2, H, T + # abd: CO, H2, H, T self.abd = {} for s in ["CO", "HCO+", "H2", "H", "e"]: - self.abd[s] = np.swapaxes( self.ytdata[s].to_ndarray(), 0, 2 ) + self.abd[s] = np.swapaxes(self.ytdata[s].to_ndarray(), 0, 2) self.abd["T"] = np.swapaxes( - (self.ytdata["temperature_chem"]/u.K).in_cgs().to_ndarray() , 0, 2 ) + (self.ytdata["temperature_chem"] / u.K).in_cgs().to_ndarray(), 0, 2 + ) self.abd["Tinit"] = np.swapaxes( - (self.ytdata["temperature_init"]/u.K).in_cgs().to_ndarray() , 0, 2 ) + (self.ytdata["temperature_init"] / u.K).in_cgs().to_ndarray(), 0, 2 + ) return + class CoreData: - def __init__(self, fn, left_edge=None, dims=None, n0=1.0e3, T0=12.): + def __init__(self, fn, left_edge=None, dims=None, n0=1.0e3, T0=12.0): self.n0 = n0 self.T0 = T0 - self.mu = 1.4 #mean molecular weight including He - self.cs0_kms = 0.19 * (T0/10.)**0.5 - L0 = 0.87 * const.pc * (n0/1.0e3)**(-0.5) * (T0/10.)**0.5 - L0_pc = L0/const.pc + self.mu = 1.4 # mean molecular weight including He + self.cs0_kms = 0.19 * (T0 / 10.0) ** 0.5 + L0 = 0.87 * const.pc * (n0 / 1.0e3) ** (-0.5) * (T0 / 10.0) ** 0.5 + L0_pc = L0 / const.pc M0 = L0**3 * n0 * const.mh * self.mu - unit_base={"length_unit": (L0_pc,"pc"), - "time_unit": (L0_pc/self.cs0_kms,"s*pc/km"), - "mass_unit": (M0,"g"), - "velocity_unit": (self.cs0_kms,"km/s")} - #load data - self.ytds = yt.load(fn, units_override = unit_base) + unit_base = { + "length_unit": (L0_pc, "pc"), + "time_unit": (L0_pc / self.cs0_kms, "s*pc/km"), + "mass_unit": (M0, "g"), + "velocity_unit": (self.cs0_kms, "km/s"), + } + # load data + self.ytds = yt.load(fn, units_override=unit_base) if left_edge is None: self.left_edge = self.ytds.domain_left_edge else: @@ -84,20 +93,21 @@ def __init__(self, fn, left_edge=None, dims=None, n0=1.0e3, T0=12.): self.dims = self.ytds.domain_dimensions else: self.dims = dims - self.ytdata = self.ytds.covering_grid(level=0, left_edge=self.left_edge, - dims=self.dims) + self.ytdata = self.ytds.covering_grid( + level=0, left_edge=self.left_edge, dims=self.dims + ) self.shape_xyz = self.dims self.nx = self.shape_xyz[0] self.ny = self.shape_xyz[1] self.nz = self.shape_xyz[2] self.dx = (self.ytdata["dx"][0, 0, 0].in_cgs()).value rho = self.ytdata["density"].in_cgs().to_ndarray() - nH0 = rho/const.mh/self.mu + nH0 = rho / const.mh / self.mu self.nH = np.swapaxes(nH0, 0, 2) - #velocity in km/s - velx0 = (self.ytdata["velocity_x"] / (u.km/u.s)).in_cgs().to_ndarray() - vely0 = (self.ytdata["velocity_y"] / (u.km/u.s)).in_cgs().to_ndarray() - velz0 = (self.ytdata["velocity_z"] / (u.km/u.s)).in_cgs().to_ndarray() + # velocity in km/s + velx0 = (self.ytdata["velocity_x"] / (u.km / u.s)).in_cgs().to_ndarray() + vely0 = (self.ytdata["velocity_y"] / (u.km / u.s)).in_cgs().to_ndarray() + velz0 = (self.ytdata["velocity_z"] / (u.km / u.s)).in_cgs().to_ndarray() vel = np.zeros([self.nx, self.ny, self.nz, 3]) for ix in range(self.nx): for iy in range(self.ny): @@ -106,35 +116,34 @@ def __init__(self, fn, left_edge=None, dims=None, n0=1.0e3, T0=12.): vel[ix, iy, iz, 1] = vely0[ix, iy, iz] vel[ix, iy, iz, 2] = velz0[ix, iy, iz] self.vel_kms = np.swapaxes(vel, 0, 2) - #abd: temprature, H2, NH3 + # abd: temprature, H2, NH3 self.abd = {} self.abd["T"] = np.zeros(self.nH.shape) + self.T0 self.abd["H2"] = np.zeros(self.nH.shape) + 0.5 xNH3 = np.zeros(self.nH.shape) n_th = 6.0e3 - xNH3[self.nH > n_th] = 10**(-8.5) + xNH3[self.nH > n_th] = 10 ** (-8.5) self.abd["pNH3"] = xNH3 return + class Data: """Read from simulation output, and functions to write RADMC-3D input files, run radmc3d to produce outputs.""" - def __init__(self, fn, fmrt="TIGRESS", left_edge=None, dims=None, - **keys): + + def __init__(self, fn, fmrt="TIGRESS", left_edge=None, dims=None, **keys): """fmrt: input format, ["TIGRESS", "core"]""" self.fn = fn self.fmrt = fmrt if self.fmrt == "TIGRESS": - self.data = YTData(self.fn, left_edge=left_edge, dims=dims, - **keys) + self.data = YTData(self.fn, left_edge=left_edge, dims=dims, **keys) for s in ["CO", "HCO+", "H2", "H", "e"]: - indx = self.data.abd[s] < 0. - self.data.abd[s][indx] = 0. + indx = self.data.abd[s] < 0.0 + self.data.abd[s][indx] = 0.0 elif self.fmrt == "core": - self.data = CoreData(self.fn, left_edge=left_edge, dims=dims, - **keys) + self.data = CoreData(self.fn, left_edge=left_edge, dims=dims, **keys) else: - raise RuntimeError("file format not recogonized: "+ self.fmrt) + raise RuntimeError("file format not recogonized: " + self.fmrt) self.dx = self.data.dx self.nx = self.data.nx self.ny = self.data.ny @@ -145,6 +154,7 @@ def __init__(self, fn, fmrt="TIGRESS", left_edge=None, dims=None, self.linemodes["LVG"] = "3" self.linemodes["Thin"] = "4" return + def get_NH2(self, char_axis): """return in [nx, ny, nz] shape""" if char_axis == "x": @@ -157,41 +167,60 @@ def get_NH2(self, char_axis): raise RuntimeError("char_axis has to be {x, y, z}") nH = np.swapaxes(self.data.nH, 0, 2) xH2 = np.swapaxes(self.data.abd["H2"], 0, 2) - return (self.dx * nH*xH2).sum(axis=nax) + return (self.dx * nH * xH2).sum(axis=nax) - def prepare(self, work_dir, image_dir, clean_outfile=True, clean_infile=True, species_coll=None, - moldata_dir="/tigress/munan/chemistry/scripts/radmc3d/prob/moldata/", - line_mode="LVG", output_binary=True, - escprob=None, escprob_max_pc=100., gas_velocity=None, micro_turbulence=None, - gas_temperature="data", gas_temperature_max=200., - temperature_bg=2.73, - gas_temperature_cold=None, slow_lvg_conv=None, - species_main="co", minimum_abundance=None, ortho2paraH2=3., - iline=1, widthkms=1, vkms=0, linenlam=40, isloadlambda=False, - external_tbg=None): + def prepare( + self, + work_dir, + image_dir, + clean_outfile=True, + clean_infile=True, + species_coll=None, + moldata_dir="/tigress/munan/chemistry/scripts/radmc3d/prob/moldata/", + line_mode="LVG", + output_binary=True, + escprob=None, + escprob_max_pc=100.0, + gas_velocity=None, + micro_turbulence=None, + gas_temperature="data", + gas_temperature_max=200.0, + temperature_bg=2.73, + gas_temperature_cold=None, + slow_lvg_conv=None, + species_main="co", + minimum_abundance=None, + ortho2paraH2=3.0, + iline=1, + widthkms=1, + vkms=0, + linenlam=40, + isloadlambda=False, + external_tbg=None, + ): """Prepare for input paramters, create simple input files. - This can be done multiple times to the Data class (work_dir). - moldata_dir: store molecular line data for co. - line_mode: LVG, LTE, Thin - escprob, gas_veolcity, micro_turbulence, gas_temperature: - {"data" (in Data class), None (not specified), float_number (constant), np.array([nx,ny,nz]) } - gas_temperature_max: cap the gas temperature at some maximum value. Default 200. Kelvin. - temperature_bg: back ground temperature. Default 2.73K for cmb. - gas_temperature_cold: upper liit of temperature of cold gas. If not - None, set species abundance in any hotter gas to be zero. - escprob_max_pc: cap the escape probability maximu, value in pc. Default 100 pc. - gas_temperature_max, gas_temperature_cold and escprob_max_pc - are only used when the choice 'data' is specified. - slow_lvg_conv: if sepecified, add slow lvg and the convergence - criteria in radmc3d.inp""" - #create working directory if it does not exist + This can be done multiple times to the Data class (work_dir). + moldata_dir: store molecular line data for co. + line_mode: LVG, LTE, Thin + escprob, gas_veolcity, micro_turbulence, gas_temperature: + {"data" (in Data class), None (not specified), float_number (constant), np.array([nx,ny,nz]) } + gas_temperature_max: cap the gas temperature at some maximum value. Default 200. Kelvin. + temperature_bg: back ground temperature. Default 2.73K for cmb. + gas_temperature_cold: upper liit of temperature of cold gas. If not + None, set species abundance in any hotter gas to be zero. + escprob_max_pc: cap the escape probability maximu, value in pc. Default 100 pc. + gas_temperature_max, gas_temperature_cold and escprob_max_pc + are only used when the choice 'data' is specified. + slow_lvg_conv: if sepecified, add slow lvg and the convergence + criteria in radmc3d.inp""" + # create working directory if it does not exist self.work_dir = work_dir self.image_dir = image_dir self.species_coll = species_coll self.species_main = species_main self.minimum_abundance = minimum_abundance - if species_coll != None: - self.ncoll = len(species_coll) #number of collisional species + if species_coll is not None: + self.ncoll = len(species_coll) # number of collisional species else: self.ncoll = 0 if not os.path.exists(work_dir): @@ -202,43 +231,53 @@ def prepare(self, work_dir, image_dir, clean_outfile=True, clean_infile=True, sp os.makedirs(image_dir) else: print("{} already exists. Using the directory.".format(image_dir)) - #clean up directory + # clean up directory if clean_outfile: print("Cleaning radmc3d output files...") cmd = "rm " + work_dir + "*\.*out " + work_dir + "*\.*dat" - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(p.stdout.read()) if clean_infile: print("Cleaning radmc3d input files...") cmd = "rm " + work_dir + "*\.*inp " - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(p.stdout.read()) - #copy molecular line data into directory + # copy molecular line data into directory print("Copy molecular data file into directory..") - fn_mol = "molecule_"+self.species_main+".inp" + fn_mol = "molecule_" + self.species_main + ".inp" cmd = "cp " + moldata_dir + fn_mol + " " + work_dir - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(p.stdout.read()) - #parameters for the line module + # parameters for the line module self.iline = iline self.widthkms = widthkms self.vkms = vkms self.linenlam = linenlam self.isloadlambda = isloadlambda self.external_tbg = external_tbg - #calculate wavelength points + # calculate wavelength points if self.isloadlambda: line_nu0 = read_line_freq(moldata_dir + fn_mol) - self.wavelength_micron = get_wavelength_micron(line_nu0, iline, widthkms, - vkms, linenlam) - #external radiation requires loadlambda to be turned on + self.wavelength_micron = get_wavelength_micron( + line_nu0, iline, widthkms, vkms, linenlam + ) + # external radiation requires loadlambda to be turned on if external_tbg is not None: if not isloadlambda: - raise RuntimeError("ERORR:" - + " external radiation option requires loadlambda to be turned on.") - external_intensity = const.bb_inu(temp_k=external_tbg, - lam_cm=self.wavelength_micron/1.0e4, units=False) - print("Creating external_source.inp..."); + raise RuntimeError( + "ERORR:" + + " external radiation option requires loadlambda to be turned on." + ) + external_intensity = const.bb_inu( + temp_k=external_tbg, lam_cm=self.wavelength_micron / 1.0e4, units=False + ) + print("Creating external_source.inp...") fw = open(work_dir + "external_source.inp", "w") fw.write("2\n{:d}".format(self.linenlam)) for li in self.wavelength_micron: @@ -246,36 +285,50 @@ def prepare(self, work_dir, image_dir, clean_outfile=True, clean_infile=True, sp for ii in external_intensity: fw.write("\n{:.6e}".format(ii)) fw.close() - #create camera_wavelength_micron.inp if loadlambda is used + # create camera_wavelength_micron.inp if loadlambda is used if self.isloadlambda: - print("Creating camera_wavelength_micron.inp required by loadlambda option..."); + print( + "Creating camera_wavelength_micron.inp required by loadlambda option..." + ) fw = open(work_dir + "camera_wavelength_micron.inp", "w") fw.write("{:d}".format(self.linenlam)) for li in self.wavelength_micron: fw.write("\n{:f}".format(li)) fw.close() - #copy to wavelength_micron.inp, used to check background radiation - cmd = "cp "+ work_dir + "camera_wavelength_micron.inp " + ( - work_dir + "wavelength_micron.inp") - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + # copy to wavelength_micron.inp, used to check background radiation + cmd = ( + "cp " + + work_dir + + "camera_wavelength_micron.inp " + + (work_dir + "wavelength_micron.inp") + ) + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(p.stdout.read()) else: - #create wavelength_micron.inp necessary for file + # create wavelength_micron.inp necessary for file print("Creating wavelength_micron.inp required by radmc3d..") fw = open(work_dir + "wavelength_micron.inp", "w") fw.write("2\n0.1\n100\n\n") fw.close() - #create lines.inp + # create lines.inp print("Writing lines.inp...") fl = open(work_dir + "lines.inp", "w") - fl.write("2\n1\n"+self.species_main+" leiden 0 0 "+str(self.ncoll) + "\n") - if self.species_coll != None: + fl.write( + "2\n1\n" + + self.species_main + + " leiden 0 0 " + + str(self.ncoll) + + "\n" + ) + if self.species_coll is not None: for s in species_coll: fl.write(s) fl.write("\n") fl.write("\n") fl.close() - #create radmc3d.inp + # create radmc3d.inp print("Writing radmc3d.inp...") fr = open(work_dir + "radmc3d.inp", "w") if output_binary: @@ -287,12 +340,16 @@ def prepare(self, work_dir, image_dir, clean_outfile=True, clean_infile=True, sp else: fr.write("lines_mode = " + self.linemodes[line_mode]) fr.write("\n") - if slow_lvg_conv != None: - fr.write("lines_nonlte_convcrit = " + str(slow_lvg_conv) + \ - "\n" + "lines_slowlvg_as_alternative=1\n") + if slow_lvg_conv is not None: + fr.write( + "lines_nonlte_convcrit = " + + str(slow_lvg_conv) + + "\n" + + "lines_slowlvg_as_alternative=1\n" + ) fr.write("\n\n") fr.close() - #store input parameters + # store input parameters self.line_mode = line_mode self.output_binary = output_binary self.escprob = escprob @@ -304,180 +361,218 @@ def prepare(self, work_dir, image_dir, clean_outfile=True, clean_infile=True, sp self.temperature_bg = temperature_bg self.gas_temperature_cold = gas_temperature_cold self.ortho2paraH2 = ortho2paraH2 - #give summary of parameters. + # give summary of parameters. print("\nSummary, radmc input parameters:") - print("species:{}, minimum abundance:{}".format(self.species_main, - self.minimum_abundance)) + print( + "species:{}, minimum abundance:{}".format( + self.species_main, self.minimum_abundance + ) + ) print("Working directory: {}".format(self.work_dir)) - print("{} collisional species: {}, oH2/pH2={}".format( - self.ncoll, self.species_coll, self.ortho2paraH2)) - #print("Note that the collisional speceis in molecule_co.inp is pH2 and oH2!") - print("line_mode={}, output_binary={}".format(self.line_mode, self.output_binary)) - print("escprob={}, escprob_max_pc={}, gas_velocity={}, micro_turbulence={}".format( - self.escprob, self.escprob_max_pc, self.gas_velocity, - self.micro_turbulence)) - print("gas_temperature={}, gas_temperature_max={}, gas_temperature_cold={}".format( - self.gas_temperature, self.gas_temperature_max, - self.gas_temperature_cold)) - print("iline={}, widthkms={}, vkms={}, linenlam={}, isloadlambda={}".format( - iline, widthkms, vkms, linenlam, isloadlambda)) + print( + "{} collisional species: {}, oH2/pH2={}".format( + self.ncoll, self.species_coll, self.ortho2paraH2 + ) + ) + # print("Note that the collisional speceis in molecule_co.inp is pH2 and oH2!") + print( + "line_mode={}, output_binary={}".format(self.line_mode, self.output_binary) + ) + print( + "escprob={}, escprob_max_pc={}, gas_velocity={}, micro_turbulence={}".format( + self.escprob, + self.escprob_max_pc, + self.gas_velocity, + self.micro_turbulence, + ) + ) + print( + "gas_temperature={}, gas_temperature_max={}, gas_temperature_cold={}".format( + self.gas_temperature, + self.gas_temperature_max, + self.gas_temperature_cold, + ) + ) + print( + "iline={}, widthkms={}, vkms={}, linenlam={}, isloadlambda={}".format( + iline, widthkms, vkms, linenlam, isloadlambda + ) + ) return def write_inputs(self): """Write main radmc3d input files *.binp""" print("\nWriting main input files for radmc3d...") - #amr_grid.binp - xmin = - self.dx * self.nx/2 - ymin = - self.dx * self.ny/2 - zmin = - self.dx * self.nz/2 - self._write_grid_binp(self.work_dir, self.nx, self.ny, self.nz, self.dx, xmin, ymin, zmin) - #numberdens_species.binp - fn_ns = self.work_dir + "numberdens_"+self.species_main+".binp" + # amr_grid.binp + xmin = -self.dx * self.nx / 2 + ymin = -self.dx * self.ny / 2 + zmin = -self.dx * self.nz / 2 + self._write_grid_binp( + self.work_dir, self.nx, self.ny, self.nz, self.dx, xmin, ymin, zmin + ) + # numberdens_species.binp + fn_ns = self.work_dir + "numberdens_" + self.species_main + ".binp" if self.species_main == "co": - ns = np.swapaxes(self.data.abd["CO"]*self.data.nH, 0, 2) + ns = np.swapaxes(self.data.abd["CO"] * self.data.nH, 0, 2) elif self.species_main == "hcop": - ns = np.swapaxes(self.data.abd["HCO+"]*self.data.nH, 0, 2) + ns = np.swapaxes(self.data.abd["HCO+"] * self.data.nH, 0, 2) elif self.species_main == "pnh3": - ns = np.swapaxes(self.data.abd["pNH3"]*self.data.nH, 0, 2) + ns = np.swapaxes(self.data.abd["pNH3"] * self.data.nH, 0, 2) elif self.species_main == "HI": - ns = np.swapaxes(self.data.abd["H"]*self.data.nH, 0, 2) + ns = np.swapaxes(self.data.abd["H"] * self.data.nH, 0, 2) else: - raise RuntimeError("species {} not implemented".format( - self.species_main)) + raise RuntimeError("species {} not implemented".format(self.species_main)) if self.gas_temperature_cold is not None: Tg1 = np.swapaxes(self.data.abd["Tinit"], 0, 2) indx_Thot = np.where(Tg1 > self.gas_temperature_cold) - ns[indx_Thot] = 0. + ns[indx_Thot] = 0.0 if self.minimum_abundance is not None: indx_empty = np.where(ns < self.minimum_abundance) - ns[indx_empty] = 0. + ns[indx_empty] = 0.0 self._write_binp(ns, self.nx, self.ny, self.nz, fn_ns) - #numberdens_***.binp for collisional species - #o-H2 and p-H2 + # numberdens_***.binp for collisional species + # o-H2 and p-H2 fo = float(self.ortho2paraH2) / float(1.0 + self.ortho2paraH2) - fp = 1. - fo - if self.species_coll != None: + fp = 1.0 - fo + if self.species_coll is not None: for si in range(self.ncoll): s = self.species_coll[si] if s == "oH2": fn_oh2 = self.work_dir + "numberdens_oH2.binp" - h20 = np.swapaxes(self.data.abd["H2"]*self.data.nH, 0, 2) + h20 = np.swapaxes(self.data.abd["H2"] * self.data.nH, 0, 2) oh20 = h20 * fo self._write_binp(oh20, self.nx, self.ny, self.nz, fn_oh2) elif s == "pH2": fn_ph2 = self.work_dir + "numberdens_pH2.binp" - h20 = np.swapaxes(self.data.abd["H2"]*self.data.nH, 0, 2) + h20 = np.swapaxes(self.data.abd["H2"] * self.data.nH, 0, 2) ph20 = h20 * fp self._write_binp(ph20, self.nx, self.ny, self.nz, fn_ph2) else: fn_s = self.work_dir + "numberdens_" + s + ".binp" - s0 = np.swapaxes(self.data.abd[s]*self.data.nH, 0, 2) + s0 = np.swapaxes(self.data.abd[s] * self.data.nH, 0, 2) self._write_binp(s0, self.nx, self.ny, self.nz, fn_s) - #co_coll = ["pH2", "oH2"] - #print("Warning: species {} used as {}".format(s, co_coll[si])) - #gas_temperature.binp + # co_coll = ["pH2", "oH2"] + # print("Warning: species {} used as {}".format(s, co_coll[si])) + # gas_temperature.binp fn_Tg = self.work_dir + "gas_temperature.binp" if self.gas_temperature == "data": Tg0 = np.swapaxes(self.data.abd["T"], 0, 2) indx_Tghigh = np.where(Tg0 > self.gas_temperature_max) Tg0[indx_Tghigh] = self.gas_temperature_max - elif type(self.gas_temperature) == float: + elif isinstance(self.gas_temperature, float): Tg0 = np.zeros([self.nx, self.ny, self.nz]) + self.gas_temperature - elif type(self.gas_temperature) == np.ndarray: + elif isinstance(self.gas_temperature, np.ndarray): if self.gas_temperature.shape == self.shape_xyz: Tg0 = self.gas_temperature else: raise RuntimeError( "gas_temperature array has shape {} != grid shape {}".format( - self.gas_temperature.shape, self.shape_xyz)) + self.gas_temperature.shape, self.shape_xyz + ) + ) else: - raise RuntimeError("gas_temperature: {}, not recognized.".format( - self.gas_temperature)) + raise RuntimeError( + "gas_temperature: {}, not recognized.".format(self.gas_temperature) + ) self._write_binp(Tg0, self.nx, self.ny, self.nz, fn_Tg) - #escprob_lengthscale.binp + # escprob_lengthscale.binp fn_Lesc = self.work_dir + "escprob_lengthscale.binp" if self.escprob == "data": - grad = np.gradient(np.swapaxes(self.data.nH, 0, 2), self.dx) - Lesc0 = 1. / np.sqrt( grad[0]**2 + grad[1]**2 + grad[2]**2) + grad = np.gradient(np.swapaxes(self.data.nH, 0, 2), self.dx) + Lesc0 = 1.0 / np.sqrt(grad[0] ** 2 + grad[1] ** 2 + grad[2] ** 2) indx = Lesc0 > self.escprob_max_pc * const.pc Lesc0[indx] = self.escprob_max_pc * const.pc - elif self.escprob == None: + elif self.escprob is None: print("No escprob_lengthscale.binp, Lesc = inf.") elif self.escprob == "jeans": - #use jeans length with temperature cap at 40K + # use jeans length with temperature cap at 40K nH0 = self.data.nH - mu = const.mh*1.4 / (1 - self.data.abd["H2"] + self.data.abd["e"] + 0.1) + mu = const.mh * 1.4 / (1 - self.data.abd["H2"] + self.data.abd["e"] + 0.1) Tcs = np.copy(self.data.abd["T"]) - indxTcs = Tcs > 40. - Tcs[indxTcs] = 40. - cs = np.sqrt( const.kb * Tcs / mu ) - LJ = cs * np.sqrt(math.pi / (const.g * nH0 * mu) ) + indxTcs = Tcs > 40.0 + Tcs[indxTcs] = 40.0 + cs = np.sqrt(const.kb * Tcs / mu) + LJ = cs * np.sqrt(math.pi / (const.g * nH0 * mu)) Lesc0 = np.swapaxes(LJ, 0, 2) - elif type(self.escprob) == float: + elif isinstance(self.escprob, float): Lesc0 = np.zeros(self.shape_xyz) + self.escprob - elif type(self.escprob) == np.ndarray: + elif isinstance(self.escprob, np.ndarray): if self.escprob.shape == self.shape_xyz: Lesc0 = self.escprob else: raise RuntimeError( "escprob array has shape {} != grid shape {}".format( - self.escprob.shape, self.shape_xyz)) + self.escprob.shape, self.shape_xyz + ) + ) else: - raise RuntimeError("escprob: {}, not recognized.".format( - self.escprob)) - if self.escprob != None: + raise RuntimeError("escprob: {}, not recognized.".format(self.escprob)) + if self.escprob is not None: self._write_binp(Lesc0, self.nx, self.ny, self.nz, fn_Lesc) - #micro_turbulence.binp + # micro_turbulence.binp fn_mt = self.work_dir + "microturbulence.binp" if self.micro_turbulence == "data": - mt0 = np.zeros(self.shape_xyz) + const.km * np.sqrt(self.dx/const.pc) - elif self.micro_turbulence == None: + mt0 = np.zeros(self.shape_xyz) + const.km * np.sqrt(self.dx / const.pc) + elif self.micro_turbulence is None: print("No micro_turbulence.binp, vturb = 0.") - elif type(self.micro_turbulence) == float: + elif isinstance(self.micro_turbulence, float): mt0 = np.zeros(self.shape_xyz) + self.micro_turbulence - elif type(self.micro_turbulence) == np.ndarray: + elif isinstance(self.micro_turbulence, np.ndarray): if self.micro_turbulence.shape == self.shape_xyz: mt0 = self.micro_turbulence else: raise RuntimeError( "micro_turbulence array has shape {} != grid shape {}".format( - self.micro_turbulence.shape, self.shape_xyz)) + self.micro_turbulence.shape, self.shape_xyz + ) + ) else: - raise RuntimeError("micro_turbulence: {}, not recognized.".format( - self.micro_turbulence)) - if self.micro_turbulence != None: + raise RuntimeError( + "micro_turbulence: {}, not recognized.".format(self.micro_turbulence) + ) + if self.micro_turbulence is not None: self._write_binp(mt0, self.nx, self.ny, self.nz, fn_mt) - #gas_velocity.binp + # gas_velocity.binp vel_shape = (self.nx, self.ny, self.nz, 3) - if self.gas_velocity == None: + if self.gas_velocity is None: print("No gas_velocity.binp, vgas = 0.") elif self.gas_velocity == "data": - fn_vel = self.work_dir + "gas_velocity.binp" + fn_vel = self.work_dir + "gas_velocity.binp" vel = self.data.vel_kms vel0 = np.swapaxes(vel, 0, 2) * 1.0e5 - elif type(self.gas_velocity) == float: + elif isinstance(self.gas_velocity, float): vel0 = np.zeros(vel_shape) + self.gas_velocity - elif type(self.gas_velocity) == np.ndarray: + elif isinstance(self.gas_velocity, np.ndarray): if self.gas_velocity.shape == vel_shape: vel0 = self.gas_velocity else: raise RuntimeError( "gas_velocity array has shape {} != required shape {}".format( - self.gas_velocity.shape, vel_shape)) + self.gas_velocity.shape, vel_shape + ) + ) else: - raise RuntimeError("gas_velocity: {}, not recognized.".format( - self.gas_velocity)) - if self.gas_velocity != None: + raise RuntimeError( + "gas_velocity: {}, not recognized.".format(self.gas_velocity) + ) + if self.gas_velocity is not None: self._write_binp(vel0, self.nx, self.ny, self.nz, fn_vel, ndim=3) return - def run(self, image_new_name=None, levelpop_new_name=None, tausurf=False, - cmd_tail="", incl=0, phi=0): + def run( + self, + image_new_name=None, + levelpop_new_name=None, + tausurf=False, + cmd_tail="", + incl=0, + phi=0, + ): """Run radmc3d in working directory, move image.dat to a new file if image_new_name is specified.""" - if levelpop_new_name != None: + if levelpop_new_name is not None: writepop = "writepop" else: writepop = "" @@ -488,38 +583,74 @@ def run(self, image_new_name=None, levelpop_new_name=None, tausurf=False, if self.isloadlambda: cmd_lambda = " loadlambda" else: - cmd_lambda = " iline {:d} widthkms ".format(self.iline) + str(self.widthkms) + \ - " vkms " + str(self.vkms) + " linenlam " + str(self.linenlam) - cmd = "cd " + self.work_dir + \ - ";radmc3d " + cmd_image + " incl "+str(incl)+" phi "+str(phi)+" npix " + str(self.nx) + \ - " sizepc " + str(self.nx*self.dx/const.pc) + cmd_lambda + " doppcatch " + \ - writepop + cmd_tail + cmd_lambda = ( + " iline {:d} widthkms ".format(self.iline) + + str(self.widthkms) + + " vkms " + + str(self.vkms) + + " linenlam " + + str(self.linenlam) + ) + cmd = ( + "cd " + + self.work_dir + + ";radmc3d " + + cmd_image + + " incl " + + str(incl) + + " phi " + + str(phi) + + " npix " + + str(self.nx) + + " sizepc " + + str(self.nx * self.dx / const.pc) + + cmd_lambda + + " doppcatch " + + writepop + + cmd_tail + ) if self.output_binary: img_ext = "bout" else: img_ext = "out" pop_ext = "dat" print(cmd) - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) - output = (p.stdout.read()).decode('ascii') + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) + output = (p.stdout.read()).decode("ascii") self.output = output fo = open(self.work_dir + "output.txt", "w") fo.write(output) fo.close() self.image_new_name = image_new_name - if image_new_name != None: - cmd = "mv " + self.work_dir+ "image." + img_ext + " " + image_new_name - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + if image_new_name is not None: + cmd = "mv " + self.work_dir + "image." + img_ext + " " + image_new_name + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(cmd) print(p.stdout.read()) - if levelpop_new_name != None: - cmd = "mv " + self.work_dir+ "levelpop_"+self.species_main+"." + pop_ext + " " + levelpop_new_name - p = Popen(cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True) + if levelpop_new_name is not None: + cmd = ( + "mv " + + self.work_dir + + "levelpop_" + + self.species_main + + "." + + pop_ext + + " " + + levelpop_new_name + ) + p = Popen( + cmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True + ) print(cmd) print(p.stdout.read()) return + def readimage(self): - #go to run directory + # go to run directory img = Images(self.image_new_name) img.radmcData = copy.deepcopy(self) return img @@ -527,23 +658,24 @@ def readimage(self): def _write_grid_binp(self, dir_out, nx, ny, nz, dx, xmin, ymin, zmin): """dx=dy=dz, cgs units, dir_out/""" fn_grid = dir_out + "amr_grid.binp" - xi_arr = np.linspace(xmin, nx*dx+xmin, nx+1) - yi_arr = np.linspace(ymin, ny*dx+ymin, ny+1) - zi_arr = np.linspace(zmin, nz*dx+zmin, nz+1) - amr_header=(1, 0, 0, 0, 1, 1, 1, nx, ny, nz) + xi_arr = np.linspace(xmin, nx * dx + xmin, nx + 1) + yi_arr = np.linspace(ymin, ny * dx + ymin, ny + 1) + zi_arr = np.linspace(zmin, nz * dx + zmin, nz + 1) + amr_header = (1, 0, 0, 0, 1, 1, 1, nx, ny, nz) fout = open(fn_grid, "wb") - amr_header_b = struct.pack('q'*len(amr_header), *amr_header) - xi_arr_b = struct.pack('d'*len(xi_arr), *xi_arr) - yi_arr_b = struct.pack('d'*len(yi_arr), *yi_arr) - zi_arr_b = struct.pack('d'*len(zi_arr), *zi_arr) + amr_header_b = struct.pack("q" * len(amr_header), *amr_header) + xi_arr_b = struct.pack("d" * len(xi_arr), *xi_arr) + yi_arr_b = struct.pack("d" * len(yi_arr), *yi_arr) + zi_arr_b = struct.pack("d" * len(zi_arr), *zi_arr) fout.write(amr_header_b) fout.write(xi_arr_b) fout.write(yi_arr_b) fout.write(zi_arr_b) fout.close() return + def _write_binp(self, arr0, nx, ny, nz, fn, ndim=1): - arr = np.zeros(nx*ny*nz*ndim) + arr = np.zeros(nx * ny * nz * ndim) i = 0 if ndim > 1: for iz in range(nz): @@ -558,47 +690,49 @@ def _write_binp(self, arr0, nx, ny, nz, fn, ndim=1): for ix in range(nx): arr[i] = arr0[ix, iy, iz] i = i + 1 - header = (1, 8, nx*ny*nz) - header_b = struct.pack('q'*len(header), *header) - arr_b = struct.pack('d'*len(arr), *arr) + header = (1, 8, nx * ny * nz) + header_b = struct.pack("q" * len(header), *header) + arr_b = struct.pack("d" * len(arr), *arr) fout1 = open(fn, "wb") fout1.write(header_b) fout1.write(arr_b) fout1.close() return + class Images: """Read RADMC-3D image output. Functions include plotting image and sed.""" + def __init__(self, fn_image, lam0=2600.75763346): """Read image from RADMC3D output. lam0: wavelength at line center in microns, used to calculate velocity v. Default: 2600.757: CO J1-0 line center.""" - ext = fn_image.split('.')[-1] + ext = fn_image.split(".")[-1] if ext == "out": self.binary_output = False elif ext == "bout": self.binary_output = True else: raise RuntimeError("Extension {} not recognized.".format(ext)) - #read image output + # read image output image_all = [] if self.binary_output: - fimage = open(fn_image, 'rb') - _ = fimage.read(8) #iformat - precis = 8 #precision - real = 'd' - lines = fimage.read(8*3) - lines = struct.unpack("q"*3, lines) + fimage = open(fn_image, "rb") + _ = fimage.read(8) # iformat + precis = 8 # precision + real = "d" + lines = fimage.read(8 * 3) + lines = struct.unpack("q" * 3, lines) nximg = lines[0] nyimg = lines[1] nlam = lines[2] lines = fimage.read(precis * 2) - lines = struct.unpack(real*2, lines) + lines = struct.unpack(real * 2, lines) pixsize_x = lines[0] pixsize_y = lines[1] lines = fimage.read(precis * nlam) - lines = struct.unpack(real*nlam, lines) + lines = struct.unpack(real * nlam, lines) lam_arr = np.array(lines) for i in range(nlam): image = np.zeros((nximg, nyimg)) @@ -607,38 +741,38 @@ def __init__(self, fn_image, lam0=2600.75763346): fi = fimage.read(precis) fi = struct.unpack(real, fi)[0] image[ix, iy] = fi - #image1 = np.swapaxes(image, 0, 1) + # image1 = np.swapaxes(image, 0, 1) image_all.append(image) fimage.close() else: - fimage = open(fn_image, 'r') - _ = fimage.readline() #format 1 - ldim = fimage.readline().split() + fimage = open(fn_image, "r") + _ = fimage.readline() # format 1 + ldim = fimage.readline().split() nximg = int(ldim[0]) nyimg = int(ldim[0]) - #print nx, ny + # print nx, ny nlam = int(fimage.readline()) - #print nlam - lpix = fimage.readline().split() + # print nlam + lpix = fimage.readline().split() pixsize_x = float(lpix[0]) pixsize_y = float(lpix[1]) - #print pixsize_x/const.pc, pixsize_y/const.pc + # print pixsize_x/const.pc, pixsize_y/const.pc lam_arr = np.zeros(nlam) for i in range(nlam): - llam = fimage.readline().split() + llam = fimage.readline().split() lam_arr[i] = float(llam[0]) - #print lam_arr + # print lam_arr for i in range(nlam): - _ = fimage.readline() #empty line + _ = fimage.readline() # empty line image = np.zeros((nximg, nyimg)) for iy in range(nyimg): for ix in range(nximg): fi = float(fimage.readline()) image[ix, iy] = fi - #image1 = np.swapaxes(image, 0, 1) + # image1 = np.swapaxes(image, 0, 1) image_all.append(image) fimage.close() - #assign to class + # assign to class self.fn = fn_image self.nlam = nlam self.lam_micron = lam_arr @@ -648,24 +782,25 @@ def __init__(self, fn_image, lam0=2600.75763346): self.dx = pixsize_x self.dy = pixsize_y self.image = image_all - self.nu = const.c / self.lam_cm #frequency in Hz - #line center wavelength and frequency + self.nu = const.c / self.lam_cm # frequency in Hz + # line center wavelength and frequency self.lam0_micron = lam0 self.lam0_cm = lam0 * 1.0e-4 self.nu0 = const.c / self.lam0_cm - #velocity relative to the line center, in km/s - self.vel_kms = - const.c * (self.nu - self.nu0)/ self.nu0 / const.km - #translate images in to Antena temperature T_A units (Kevin) + # velocity relative to the line center, in km/s + self.vel_kms = -const.c * (self.nu - self.nu0) / self.nu0 / const.km + # translate images in to Antena temperature T_A units (Kevin) self.image_TA = [] for i in range(nlam): im = self.image[i] nui = self.nu[i] - fac_TA = const.c**2 / (2. * const.kb * nui**2) + fac_TA = const.c**2 / (2.0 * const.kb * nui**2) self.image_TA.append(im * fac_TA) - #get the average image + # get the average image self.image_avg = [im.mean() for im in self.image] self.image_TA_avg = [im.mean() for im in self.image_TA] return + def clip(self, ixmin, ixmax, iymin, iymax): self.nx = ixmax - ixmin self.ny = iymax - iymin @@ -682,17 +817,23 @@ def clip(self, ixmin, ixmax, iymin, iymax): self.image_avg = [im.mean() for im in self.image] self.image_TA_avg = [im.mean() for im in self.image_TA] return - def write_fits(self, fn_fits="image.fits", - clip_range=None, comment="", axis_name1="x", axis_name2="y", - extent=None): + + def write_fits( + self, + fn_fits="image.fits", + clip_range=None, + comment="", + axis_name1="x", + axis_name2="y", + extent=None, + ): if clip_range is not None: - self.clip(clip_range[0], clip_range[1], - clip_range[2], clip_range[3]) + self.clip(clip_range[0], clip_range[1], clip_range[2], clip_range[3]) if extent is None: - xmin = -self.dx*self.nx/2./const.pc - xmax = self.dx*self.nx/2./const.pc - ymin = -self.dy*self.ny/2./const.pc - ymax = self.dy*self.ny/2./const.pc + xmin = -self.dx * self.nx / 2.0 / const.pc + xmax = self.dx * self.nx / 2.0 / const.pc + ymin = -self.dy * self.ny / 2.0 / const.pc + ymax = self.dy * self.ny / 2.0 / const.pc else: xmin = extent[0] xmax = extent[1] @@ -702,24 +843,24 @@ def write_fits(self, fn_fits="image.fits", primary_hdu = fits.PrimaryHDU() primary_hdu.header.set("author", "Munan Gong") primary_hdu.header.set("comment", comment) - primary_hdu.header.set(axis_name1.upper()+"MIN", xmin) - primary_hdu.header.set(axis_name1.upper()+"MAX", xmax) - primary_hdu.header.set("D"+axis_name1.upper(), self.dx/const.pc) - primary_hdu.header.set(axis_name2.upper()+"MIN", ymin) - primary_hdu.header.set(axis_name2.upper()+"MAX", ymax) - primary_hdu.header.set("D"+axis_name2.upper(), self.dy/const.pc) + primary_hdu.header.set(axis_name1.upper() + "MIN", xmin) + primary_hdu.header.set(axis_name1.upper() + "MAX", xmax) + primary_hdu.header.set("D" + axis_name1.upper(), self.dx / const.pc) + primary_hdu.header.set(axis_name2.upper() + "MIN", ymin) + primary_hdu.header.set(axis_name2.upper() + "MAX", ymax) + primary_hdu.header.set("D" + axis_name2.upper(), self.dy / const.pc) primary_hdu.header.set("VMIN", self.vel_kms[0]) primary_hdu.header.set("VMAX", self.vel_kms[-1]) primary_hdu.header.set("DV", self.vel_kms[1] - self.vel_kms[0]) hdu_data = fits.ImageHDU(data) hdu_data.header.set("EXTNAME", "TA") for hdui in [primary_hdu, hdu_data]: - hdui.header.set("CDELT1", self.dx/const.pc) + hdui.header.set("CDELT1", self.dx / const.pc) hdui.header.set("CTYPE1", axis_name1) hdui.header.set("CUNIT1", "pc") hdui.header.set("CRVAL1", xmin) hdui.header.set("CRPIX1", 0) - hdui.header.set("CDELT2", self.dy/const.pc) + hdui.header.set("CDELT2", self.dy / const.pc) hdui.header.set("CTYPE2", axis_name2) hdui.header.set("CUNIT2", "pc") hdui.header.set("CRVAL2", ymin) @@ -736,16 +877,19 @@ def write_fits(self, fn_fits="image.fits", def plotimage(self, ax, TA=True, ilam=0, **keys): """If TA=True, show in Antena temperature Kelvin units. Otherwise, show in orignal units of erg/s/cm^2/Hz/str""" - if TA == True: + if TA: imgi = self.image_TA[ilam] clabel = r"$T_A (\mathrm{K})$" else: imgi = self.image[ilam] clabel = r"$I_\nu (\mathrm{erg/s/cm^2/Hz/ster})$" - cax = ax.imshow(np.swapaxes(imgi, 0, 1) + 1e-100, origin="lower", - extent=[0, self.nx*self.dx/const.pc, - 0, self.ny*self.dy/const.pc], - cmap="magma", **keys) + cax = ax.imshow( + np.swapaxes(imgi, 0, 1) + 1e-100, + origin="lower", + extent=[0, self.nx * self.dx / const.pc, 0, self.ny * self.dy / const.pc], + cmap="magma", + **keys, + ) cbar = ax.figure.colorbar(cax) cbar.solids.set_edgecolor("face") cbar.ax.set_ylabel(clabel) @@ -760,66 +904,79 @@ def plotsed(self, ax, **keys): ax.set_xlabel(r"$v(\mathrm{km/s})$") ax.set_ylabel(r"$\langle T_A \rangle (\mathrm{K})$") return + def getavgWCO(self, Tdect=0.6): """Calculate the WCO averaged over the whole image. Tdect: the detection limit in Kevin. Defaut 0.6 (K).""" - #calculate the moment 0 image + # calculate the moment 0 image img = self.getimageWCO(Tdect=Tdect) return img.mean() + def getimageWCO(self, Tdect=0.6, dv=0.07): """Return image in K km/s.""" - #get more finely spaced velocity array - #spacing in velocity array + # get more finely spaced velocity array + # spacing in velocity array vel0 = self.vel_kms vel_arr = np.arange(vel0[0], vel0[-1], dv) imageWCO = np.zeros((self.nx, self.ny)) dv0 = np.diff(vel0)[0] - #if small velocity channel, interpolation + # if small velocity channel, interpolation if dv <= dv0: - #loop over each pixel + # loop over each pixel for ix in range(self.nx): for iy in range(self.ny): - TA_arr = np.array( [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + TA_arr = np.array( + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) TA_dect = np.interp(vel_arr, vel0, TA_arr) indx = TA_dect < Tdect TA_dect[indx] = 0 - #integration + # integration imageWCO[ix, iy] = TA_dect.sum() * dv - #if large velocity channel, need to calculate flux in each channel + # if large velocity channel, need to calculate flux in each channel else: - dv_fine = min(dv/10., dv0) + dv_fine = min(dv / 10.0, dv0) vel_arr_fine = np.arange(vel0[0], vel0[-1], dv_fine) for ix in range(self.nx): for iy in range(self.ny): - TA_arr = np.array( [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + TA_arr = np.array( + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) TA_dect_fine = np.interp(vel_arr_fine, vel0, TA_arr) TA_dect = np.zeros(len(vel_arr)) for i in range(len(vel_arr)): - if i < len(vel_arr)-1: + if i < len(vel_arr) - 1: indx_channel = (vel_arr_fine >= vel_arr[i]) & ( - vel_arr_fine < vel_arr[i+1]) + vel_arr_fine < vel_arr[i + 1] + ) else: indx_channel = vel_arr_fine >= vel_arr[i] - TA_dect[i] = TA_dect_fine[indx_channel].sum() *dv_fine / dv + TA_dect[i] = TA_dect_fine[indx_channel].sum() * dv_fine / dv indx = TA_dect < Tdect TA_dect[indx] = 0 - #integration + # integration imageWCO[ix, iy] = TA_dect.sum() * dv return imageWCO - def getTexc(self, Tdect=0., dv=None): + + def getTexc(self, Tdect=0.0, dv=None): Tmax_axis = self.getTpeak(Tdect=Tdect, dv=dv) - Texc = 5.5/np.log(1+5.5/(Tmax_axis)) + Texc = 5.5 / np.log(1 + 5.5 / (Tmax_axis)) return Texc - def getTpeak(self, Tdect=0., dv=None): - if dv is None: - nv = self.nlam - else: + + def getTpeak(self, Tdect=0.0, dv=None): + if dv is not None: + # nv = self.nlam + # else: vel0 = self.vel_kms dv0 = np.diff(vel0)[0] vel_arr = np.arange(vel0[0], vel0[-1], dv) - nv = len(vel_arr) + # nv = len(vel_arr) Tmax = np.zeros((self.nx, self.ny)) for ix in range(self.nx): for iy in range(self.ny): @@ -827,41 +984,62 @@ def getTpeak(self, Tdect=0., dv=None): ITA = np.array([m[ix, iy] for m in self.image_TA]) elif dv < dv0: TA_arr = np.array( - [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) ITA = np.interp(vel_arr, vel0, TA_arr) else: vel_arr = np.arange(vel0[0], vel0[-1], dv) - dv_fine = min(dv/10., dv0) + dv_fine = min(dv / 10.0, dv0) vel_arr_fine = np.arange(vel0[0], vel0[-1], dv_fine) TA_arr = np.array( - [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) TA_dect_fine = np.interp(vel_arr_fine, vel0, TA_arr) ITA = np.zeros(len(vel_arr)) for i in range(len(vel_arr)): - if i < len(vel_arr)-1: + if i < len(vel_arr) - 1: indx_channel = (vel_arr_fine >= vel_arr[i]) & ( - vel_arr_fine < vel_arr[i+1]) + vel_arr_fine < vel_arr[i + 1] + ) else: indx_channel = vel_arr_fine >= vel_arr[i] - ITA[i] = TA_dect_fine[ - indx_channel].sum() *dv_fine / dv + ITA[i] = TA_dect_fine[indx_channel].sum() * dv_fine / dv indx_nodect = ITA < Tdect - ITA[indx_nodect] = 0. + ITA[indx_nodect] = 0.0 Tmax[ix, iy] = ITA.max() return Tmax - def plotimageWCO(self, ax, Tdect=0.6, extent=None, - cax=None,orientation="vertical", cbticks=None, **keys): + + def plotimageWCO( + self, + ax, + Tdect=0.6, + extent=None, + cax=None, + orientation="vertical", + cbticks=None, + **keys, + ): img = self.getimageWCO(Tdect=Tdect) - if extent == None: - extent=[0, self.nx*self.dx/const.pc/1e3, - 0, self.ny*self.dy/const.pc/1e3] - cax1 = ax.imshow(np.swapaxes(img, 0, 1)+1e-100, origin="lower", extent=extent, - **keys) - if cbticks != None: - cbar = ax.figure.colorbar(cax1, cax=cax,ticks=cbticks, - orientation=orientation) + if extent is None: + extent = [ + 0, + self.nx * self.dx / const.pc / 1e3, + 0, + self.ny * self.dy / const.pc / 1e3, + ] + cax1 = ax.imshow( + np.swapaxes(img, 0, 1) + 1e-100, origin="lower", extent=extent, **keys + ) + if cbticks is not None: + cbar = ax.figure.colorbar( + cax1, cax=cax, ticks=cbticks, orientation=orientation + ) else: cbar = ax.figure.colorbar(cax1, cax=cax, orientation=orientation) cbar.solids.set_edgecolor("face") @@ -869,33 +1047,52 @@ def plotimageWCO(self, ax, Tdect=0.6, extent=None, ax.set_xlabel("$x(\mathrm{kpc})$", fontsize=20) ax.set_ylabel("$y(\mathrm{kpc})$", fontsize=20) return - def plotimageXCO(self, ax, imgNH2, Tdect=0.6, WCOmin=0.6, extent=None, - cbticks=None, cax=None, orientation="vertical", **keys): + + def plotimageXCO( + self, + ax, + imgNH2, + Tdect=0.6, + WCOmin=0.6, + extent=None, + cbticks=None, + cax=None, + orientation="vertical", + **keys, + ): """WCOmin: below which assign XCO=0 (non-detection.)""" img = self.getimageWCO(Tdect=Tdect) - imgXCO = imgNH2/2.0e20/img + imgXCO = imgNH2 / 2.0e20 / img indx_WCO = img < WCOmin imgXCO[indx_WCO] = 0 - if extent == None: - extent=[0, self.nx*self.dx/const.pc/1e3, - 0, self.ny*self.dy/const.pc/1e3] - cax1 = ax.imshow(np.swapaxes(imgXCO, 0, 1), origin="lower", - extent=extent, **keys) - if cbticks != None: - cbar = ax.figure.colorbar(cax1, cax=cax,ticks=cbticks, - orientation=orientation) + if extent is None: + extent = [ + 0, + self.nx * self.dx / const.pc / 1e3, + 0, + self.ny * self.dy / const.pc / 1e3, + ] + cax1 = ax.imshow( + np.swapaxes(imgXCO, 0, 1), origin="lower", extent=extent, **keys + ) + if cbticks is not None: + cbar = ax.figure.colorbar( + cax1, cax=cax, ticks=cbticks, orientation=orientation + ) else: - cbar = ax.figure.colorbar(cax1, cax=cax ,orientation=orientation) + cbar = ax.figure.colorbar(cax1, cax=cax, orientation=orientation) cbar.solids.set_edgecolor("face") cbar.set_label( r"$\mathrm{X_\mathrm{CO}}(\mathrm{2\times 10^{20}cm^{-2}K^{-1} km^{-1}s})$", - fontsize=25) + fontsize=25, + ) ax.set_xlabel("$x(\mathrm{kpc})$", fontsize=20) ax.set_ylabel("$y(\mathrm{kpc})$", fontsize=20) return - def getimage_velocities(self, Tdect=0., dv=None): + + def getimage_velocities(self, Tdect=0.0, dv=None): """img_v_avg_kms: intensty weighted mean velocity, in unit of km/s. - img_sigma_v: intensity weighted veolcity dispersion, in unit of km/s.""" + img_sigma_v: intensity weighted veolcity dispersion, in unit of km/s.""" img_v_avg_kms = np.zeros((self.nx, self.ny)) img_sigma_v_kms = np.zeros((self.nx, self.ny)) vel0 = self.vel_kms @@ -906,36 +1103,42 @@ def getimage_velocities(self, Tdect=0., dv=None): ITA = np.array([m[ix, iy] for m in self.image_TA]) vel_arr = self.vel_kms elif dv < dv0: - #if small velocity channel, interpolation + # if small velocity channel, interpolation vel_arr = np.arange(vel0[0], vel0[-1], dv) TA_arr = np.array( - [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) ITA = np.interp(vel_arr, vel0, TA_arr) else: - #if large velocity channel, - #need to calculate flux in each channel + # if large velocity channel, + # need to calculate flux in each channel vel_arr = np.arange(vel0[0], vel0[-1], dv) - dv_fine = min(dv/10., dv0) + dv_fine = min(dv / 10.0, dv0) vel_arr_fine = np.arange(vel0[0], vel0[-1], dv_fine) TA_arr = np.array( - [(self.image_TA[self.nlam-ilam-1])[ix, iy] - for ilam in np.arange(self.nlam)]) + [ + (self.image_TA[self.nlam - ilam - 1])[ix, iy] + for ilam in np.arange(self.nlam) + ] + ) TA_dect_fine = np.interp(vel_arr_fine, vel0, TA_arr) ITA = np.zeros(len(vel_arr)) for i in range(len(vel_arr)): - if i < len(vel_arr)-1: + if i < len(vel_arr) - 1: indx_channel = (vel_arr_fine >= vel_arr[i]) & ( - vel_arr_fine < vel_arr[i+1]) + vel_arr_fine < vel_arr[i + 1] + ) else: indx_channel = vel_arr_fine >= vel_arr[i] - ITA[i] = TA_dect_fine[ - indx_channel].sum() *dv_fine / dv + ITA[i] = TA_dect_fine[indx_channel].sum() * dv_fine / dv indx_nodect = ITA < Tdect - ITA[indx_nodect] = 0. - if ITA.sum() == 0.: - img_v_avg_kms[ix, iy] = 0. - img_sigma_v_kms[ix, iy] = 0. + ITA[indx_nodect] = 0.0 + if ITA.sum() == 0.0: + img_v_avg_kms[ix, iy] = 0.0 + img_sigma_v_kms[ix, iy] = 0.0 else: v_avg = np.average(vel_arr, weights=ITA) vsq_avg = np.average(vel_arr**2, weights=ITA) @@ -943,16 +1146,17 @@ def getimage_velocities(self, Tdect=0., dv=None): img_v_avg_kms[ix, iy] = v_avg img_sigma_v_kms[ix, iy] = sigma_v return img_v_avg_kms, img_sigma_v_kms + def rescale(self, factor, circular_beam=False, stamp=None): """Rescale image to resolution = orignal_resolution*factor. This is to see the effect of a larger observational beam.""" - #check whether nx and ny is dividable by factor + # check whether nx and ny is dividable by factor if (self.nx % factor) != 0 or (self.ny % factor) != 0: print("ERROR: nx or ny not dividable by factor") return - #make a copy of the class + # make a copy of the class img_copy = copy.deepcopy(self) - #rescale image + # rescale image img_copy.nx = int(self.nx / factor) img_copy.ny = int(self.ny / factor) img_copy.dx = self.dx * factor @@ -961,27 +1165,27 @@ def rescale(self, factor, circular_beam=False, stamp=None): img_copy.image_TA = [] for i in range(self.nlam): im = self.image[i] - im0 = im.reshape([img_copy.nx, factor, img_copy.ny, - factor]).mean(3).mean(1) + im0 = im.reshape([img_copy.nx, factor, img_copy.ny, factor]).mean(3).mean(1) if circular_beam: im1 = mcb.stamp_avg(im0, stamp) else: im1 = im0 img_copy.image.append(im1) nui = self.nu[i] - fac_TA = const.c**2 / (2. * const.kb * nui**2) + fac_TA = const.c**2 / (2.0 * const.kb * nui**2) img_copy.image_TA.append(im1 * fac_TA) img_copy.image_avg = [im.mean() for im in img_copy.image] img_copy.image_TA_avg = [im.mean() for im in img_copy.image_TA] return img_copy + class LevelPop: def __init__(self, fn, dims=None): """dims = [ny, nx, nz], if not given, then 1D array, assuming ASCII output""" self.fn = fn f = open(fn, "r") - _ = f.readline() #skip first line + _ = f.readline() # skip first line self.ncells = int(f.readline()) self.nlevels = int(f.readline()) self.levels = [int(x) for x in f.readline().split()] @@ -989,13 +1193,15 @@ def __init__(self, fn, dims=None): raise RuntimeError("File: " + fn + ", nlevels != len(levels)") self.levelpop = {} for ilevel in self.levels: - if dims == None: + if dims is None: self.levelpop[ilevel] = np.zeros(self.ncells) else: if dims[0] * dims[1] * dims[2] != self.ncells or len(dims) != 3: - raise RuntimeError("File: " + fn + ", dims is not consistent with ncells.") + raise RuntimeError( + "File: " + fn + ", dims is not consistent with ncells." + ) self.levelpop[ilevel] = np.zeros(dims) - if dims == None: + if dims is None: for i in range(self.ncells): li = [float(x) for x in f.readline().split()] for ii, ilevel in zip(range(self.nlevels), self.levels): @@ -1009,33 +1215,36 @@ def __init__(self, fn, dims=None): (self.levelpop[ilevel])[ix, iy, iz] = li[ii] f.close() return + def getTexc(self, Jl=1, Ju=2, Tul=5.5, gl=1, gu=3): """Return excitation temperature. - Note that in RADMC, Ju is defined to start from 1. + Note that in RADMC, Ju is defined to start from 1. """ nl = self.levelpop[Jl] nu = self.levelpop[Ju] + 1e-100 - Texc = Tul / np.log((nl/float(gl)) / (nu/float(gu))) + Texc = Tul / np.log((nl / float(gl)) / (nu / float(gu))) return Texc + def bplanck(temp, nu): """Black body planck function B_nu(T), copied from radmc source code. - 2 h nu^3 / c^2 - B_nu(T) = ------------------ [ erg / cm^2 s ster Hz ] - exp(h nu / kT) - 1 + 2 h nu^3 / c^2 + B_nu(T) = ------------------ [ erg / cm^2 s ster Hz ] + exp(h nu / kT) - 1 - ARGUMENTS: - temp [K] = Temperature - nu [Hz] = Frequency + ARGUMENTS: + temp [K] = Temperature + nu [Hz] = Frequency """ small_ = 1.0e-200 - large_ = 709.78271 - if (temp < small_): - return 0. + # large_ = 709.78271 + if temp < small_: + return 0.0 xx = 4.7989e-11 * nu / temp - ret = 1.47455e-47 * nu * nu * nu / (np.exp(xx)-1.) + small_ + ret = 1.47455e-47 * nu * nu * nu / (np.exp(xx) - 1.0) + small_ return ret + def read_line_freq(fn_mol): """Read frequencies for molecular lines. input: @@ -1046,7 +1255,7 @@ def read_line_freq(fn_mol): f = open(fn_mol) lines = f.readlines() nlevel = int(lines[5].strip()) - ntrans = int(lines[8+nlevel].strip()) + ntrans = int(lines[8 + nlevel].strip()) line_nu0 = np.zeros(ntrans) for i in np.arange(ntrans): indx = i + 10 + nlevel @@ -1054,11 +1263,15 @@ def read_line_freq(fn_mol): line_nu0[i] = float(line_list[4]) * 1.0e9 return line_nu0 + def get_wavelength_micron(line_nu0, iline, widthkms, vkms, linenlam): - nu0 = line_nu0[iline-1] + nu0 = line_nu0[iline - 1] camera_nu = np.zeros(linenlam) for inu in np.arange(linenlam): - camera_nu[inu] = nu0 * (1. - vkms*const.km/const.c - - (2.*inu/(linenlam-1.)-1.) * widthkms*const.km/const.c) + camera_nu[inu] = nu0 * ( + 1.0 + - vkms * const.km / const.c + - (2.0 * inu / (linenlam - 1.0) - 1.0) * widthkms * const.km / const.c + ) wavelength_micron = 1.0e4 * const.c / camera_nu return wavelength_micron diff --git a/astro_tigress/synthetic_HI.py b/astro_tigress/synthetic_HI.py index 8a39630..54c2679 100644 --- a/astro_tigress/synthetic_HI.py +++ b/astro_tigress/synthetic_HI.py @@ -1,104 +1,118 @@ import numpy as np import tqdm -def los_to_HI_axis_proj(dens,temp,vel,vchannel,deltas=1.,los_axis=1): +# save synthetic HI data to fits +from astropy.io import fits + + +def los_to_HI_axis_proj(dens, temp, vel, vchannel, deltas=1.0, los_axis=1): """ - inputs: - dens: number density of hydrogen in units of 1/cm^3 - temp: temperature in units of K - vel: line-of-sight velocity in units of km/s - vchannel: velocity channel in km/s - parameters: - deltas: length of line segments in units of pc - memlim: memory limit in GB - los_axis: 0 -- z, 1 -- y, 2 -- x - outputs: a dictionary - TB: the brightness temperature - tau: optical depth + inputs: + dens: number density of hydrogen in units of 1/cm^3 + temp: temperature in units of K + vel: line-of-sight velocity in units of km/s + vchannel: velocity channel in km/s + parameters: + deltas: length of line segments in units of pc + memlim: memory limit in GB + los_axis: 0 -- z, 1 -- y, 2 -- x + outputs: a dictionary + TB: the brightness temperature + tau: optical depth """ - Nv=len(vchannel) + # Nv=len(vchannel) - Tlos=temp - vlos=vel - nlos=dens + Tlos = temp + vlos = vel + nlos = dens - Tspin=Tlos + Tspin = Tlos - ds=deltas*3.085677581467192e+18 + ds = deltas * 3.085677581467192e18 - v_L=0.21394414*np.sqrt(Tlos) # in units of km/s + v_L = 0.21394414 * np.sqrt(Tlos) # in units of km/s - TB=[] - tau_v=[] + TB = [] + tau_v = [] for vch in tqdm.tqdm(vchannel): - phi_v=0.00019827867/v_L*np.exp(-(1.6651092223153954* - (vch-vlos)/v_L)**2) # time - kappa_v=2.6137475e-15*nlos/Tspin*phi_v # area/volume = 1/length - tau_los=kappa_v*ds # dimensionless + phi_v = ( + 0.00019827867 + / v_L + * np.exp(-((1.6651092223153954 * (vch - vlos) / v_L) ** 2)) + ) # time + kappa_v = 2.6137475e-15 * nlos / Tspin * phi_v # area/volume = 1/length + tau_los = kappa_v * ds # dimensionless - tau_cumul=tau_los.cumsum(axis=los_axis) + tau_cumul = tau_los.cumsum(axis=los_axis) # same unit with Tspin - TB.append(np.nansum(Tspin*(1-np.exp(-tau_los))*np.exp(-tau_cumul),axis=los_axis)) + TB.append( + np.nansum( + Tspin * (1 - np.exp(-tau_los)) * np.exp(-tau_cumul), axis=los_axis + ) + ) # dimensionless - tau_v.append(np.nansum(kappa_v*ds,axis=los_axis)) - return np.array(TB),np.array(tau_v) + tau_v.append(np.nansum(kappa_v * ds, axis=los_axis)) + return np.array(TB), np.array(tau_v) def k10h(T2): """ - input: T/100K - output: collisional excitation rate of hydrogen in c.g.s. + input: T/100K + output: collisional excitation rate of hydrogen in c.g.s. """ - k10_1=1.19e-10*T2**(0.74-0.20*np.log(T2)) - k10_2=2.24e-10*T2**0.207*np.exp(-0.876/T2) - k10=k10_1 - idx=k10_2 > k10_1 - k10[idx]=k10_2[idx] + k10_1 = 1.19e-10 * T2 ** (0.74 - 0.20 * np.log(T2)) + k10_2 = 2.24e-10 * T2**0.207 * np.exp(-0.876 / T2) + k10 = k10_1 + idx = k10_2 > k10_1 + k10[idx] = k10_2[idx] return k10 + def k10e(T2): """ - input: T/100K - output: collisional excitation rate of electrion in c.g.s. + input: T/100K + output: collisional excitation rate of electrion in c.g.s. """ - Temp=T2*1.e2 - k10=-9.607+0.5*np.log10(Temp)*np.exp(-(np.log10(Temp))**4.5/1.8e3) - return 10.**k10 + Temp = T2 * 1.0e2 + k10 = -9.607 + 0.5 * np.log10(Temp) * np.exp(-((np.log10(Temp)) ** 4.5) / 1.8e3) + return 10.0**k10 -# save synthetic HI data to fits -from astropy.io import fits def create_fits(ytds): - tMyr=ytds.current_time.to('Myr').value - le=ytds.domain_left_edge.to('pc').value - re=ytds.domain_right_edge.to('pc').value - dx=(ytds.domain_width/ytds.domain_dimensions).to('pc').value + tMyr = ytds.current_time.to("Myr").value + le = ytds.domain_left_edge.to("pc").value + re = ytds.domain_right_edge.to("pc").value + dx = (ytds.domain_width / ytds.domain_dimensions).to("pc").value hdr = fits.Header() - hdr['time']=float(tMyr) - hdr['xmin']=(le[0],'pc') - hdr['xmax']=(re[0],'pc') - hdr['ymin']=(le[1],'pc') - hdr['ymax']=(re[1],'pc') - hdr['zmin']=(le[2],'pc') - hdr['zmax']=(re[2],'pc') - hdr['dx']=(dx[0],'pc') - hdr['dy']=(dx[1],'pc') - hdr['dz']=(dx[2],'pc') + hdr["time"] = float(tMyr) + hdr["xmin"] = (le[0], "pc") + hdr["xmax"] = (re[0], "pc") + hdr["ymin"] = (le[1], "pc") + hdr["ymax"] = (re[1], "pc") + hdr["zmin"] = (le[2], "pc") + hdr["zmax"] = (re[2], "pc") + hdr["dx"] = (dx[0], "pc") + hdr["dy"] = (dx[1], "pc") + hdr["dz"] = (dx[2], "pc") hdu = fits.PrimaryHDU(header=hdr) return hdu -def add_header_for_glue(hdu,hdr,axis='xyz'): - for i,ax in enumerate(axis): - hdu.header['CDELT{}'.format(i+1)]=hdr['d{}'.format(ax)] - hdu.header['CTYPE{}'.format(i+1)]=ax - hdu.header['CUNIT{}'.format(i+1)]=hdr.comments['d{}'.format(ax)] - hdu.header['CRVAL{}'.format(i+1)]=hdr['{}min'.format(ax)] - hdu.header['CRPIX{}'.format(i+1)]=hdr['{}max'.format(ax)]+hdr['{}min'.format(ax)] + +def add_header_for_glue(hdu, hdr, axis="xyz"): + for i, ax in enumerate(axis): + hdu.header["CDELT{}".format(i + 1)] = hdr["d{}".format(ax)] + hdu.header["CTYPE{}".format(i + 1)] = ax + hdu.header["CUNIT{}".format(i + 1)] = hdr.comments["d{}".format(ax)] + hdu.header["CRVAL{}".format(i + 1)] = hdr["{}min".format(ax)] + hdu.header["CRPIX{}".format(i + 1)] = ( + hdr["{}max".format(ax)] + hdr["{}min".format(ax)] + ) return -def save_to_fits(ytds,vchannel,TB,tau,axis='xyv',fitsname=None): + +def save_to_fits(ytds, vchannel, TB, tau, axis="xyv", fitsname=None): """Function warpper for fits writing Parameters @@ -119,19 +133,19 @@ def save_to_fits(ytds,vchannel,TB,tau,axis='xyv',fitsname=None): hdul = fits.HDUList() hdu = create_fits(ytds) - hdu.header['vmin']=(vchannel.min(),'km/s') - hdu.header['vmax']=(vchannel.max(),'km/s') - hdu.header['dv']=(vchannel[1]-vchannel[0],'km/s') + hdu.header["vmin"] = (vchannel.min(), "km/s") + hdu.header["vmax"] = (vchannel.max(), "km/s") + hdu.header["dv"] = (vchannel[1] - vchannel[0], "km/s") hdul.append(hdu) - for fdata,label in zip([TB,tau],['TB','tau']): - hdul.append(fits.ImageHDU(name=label,data=fdata)) + for fdata, label in zip([TB, tau], ["TB", "tau"]): + hdul.append(fits.ImageHDU(name=label, data=fdata)) - hdr=hdu.header + hdr = hdu.header for hdu in hdul: - add_header_for_glue(hdu,hdr,axis=axis) + add_header_for_glue(hdu, hdr, axis=axis) if fitsname is not None: - hdul.writeto(fitsname,overwrite=True) + hdul.writeto(fitsname, overwrite=True) print("FITS file named {} is created".format(fitsname)) return hdul diff --git a/astro_tigress/synthetic_dustpol.py b/astro_tigress/synthetic_dustpol.py index fb1f6b8..04261eb 100644 --- a/astro_tigress/synthetic_dustpol.py +++ b/astro_tigress/synthetic_dustpol.py @@ -1,43 +1,43 @@ -def calc_IQU(nH,Bx,By,Bz,deltas,los='y'): - ''' - Flat sky projection of dust polarization at 353 GHz - Parameters assumed are identical to Kim, Choi, and Flauger (2019) - Tdust = 18K - sigma_d,353 = 1.2e-26 cm^2 - p0 = 0.2 - ''' - if los=='x': +def calc_IQU(nH, Bx, By, Bz, deltas, los="y"): + """ + Flat sky projection of dust polarization at 353 GHz + Parameters assumed are identical to Kim, Choi, and Flauger (2019) + Tdust = 18K + sigma_d,353 = 1.2e-26 cm^2 + p0 = 0.2 + """ + if los == "x": Blos = Bx Bperp_x = By Bperp_y = Bz ilos = 2 - elif los=='y': + elif los == "y": Blos = By Bperp_x = -Bx Bperp_y = Bz ilos = 1 - elif los=='z': + elif los == "z": Blos = Bz Bperp_x = Bx Bperp_y = By ilos = 0 - args={'Bnu':41495.876171482356, 'sigma':1.2e-26, 'p0':0.2} - Bnu=args['Bnu'] - p0=args['p0'] - sigma=args['sigma'] + args = {"Bnu": 41495.876171482356, "sigma": 1.2e-26, "p0": 0.2} + Bnu = args["Bnu"] + p0 = args["p0"] + sigma = args["sigma"] - Bperp2=Bperp_x*Bperp_x+Bperp_y*Bperp_y - B2=Bperp2+Blos*Blos - cos2phi=(Bperp_y*Bperp_y-Bperp_x*Bperp_x)/Bperp2 - sin2phi=-Bperp_x*Bperp_y/Bperp2*2 - cosgam2=Bperp2/B2 + Bperp2 = Bperp_x * Bperp_x + Bperp_y * Bperp_y + B2 = Bperp2 + Blos * Blos + cos2phi = (Bperp_y * Bperp_y - Bperp_x * Bperp_x) / Bperp2 + sin2phi = -Bperp_x * Bperp_y / Bperp2 * 2 + cosgam2 = Bperp2 / B2 - ds=deltas*3.085677581467192e+18 - dtau=sigma*nH*ds + ds = deltas * 3.085677581467192e18 + dtau = sigma * nH * ds - I=Bnu*(1.0-p0*(cosgam2-2./3.0))*dtau - Q=p0*Bnu*cos2phi*cosgam2*dtau - U=p0*Bnu*sin2phi*cosgam2*dtau + Inu = Bnu * (1.0 - p0 * (cosgam2 - 2.0 / 3.0)) * dtau + Qnu = p0 * Bnu * cos2phi * cosgam2 * dtau + Unu = p0 * Bnu * sin2phi * cosgam2 * dtau - return I.sum(axis=ilos),Q.sum(axis=ilos),U.sum(axis=ilos) + return Inu.sum(axis=ilos), Qnu.sum(axis=ilos), Unu.sum(axis=ilos) diff --git a/astro_tigress/ytathena.py b/astro_tigress/ytathena.py index 3e08be5..dcf2afb 100644 --- a/astro_tigress/ytathena.py +++ b/astro_tigress/ytathena.py @@ -1,27 +1,22 @@ import matplotlib.pyplot as plt import yt -from yt.data_objects.level_sets.api import * +from yt.data_objects.level_sets.api import Clump import yt.units as u -from yt import add_field from yt import YTQuantity -from yt.utilities.physical_constants import \ - mh, \ - me, \ - sigma_thompson, \ - clight, \ - kboltz, \ - G +from yt.utilities.physical_constants import mh, kboltz import numpy as np import copy import math from yt.funcs import mylog import os +import const + dirpath = os.path.dirname(__file__) mylog.setLevel(50) global Zdg, GxHe, GxC, GxO, GxSi -Gzdg = 1. +Gzdg = 1.0 GxHe = 0.1 GxC = 1.6e-4 * Gzdg GxO = 3.2e-4 * Gzdg @@ -32,7 +27,7 @@ class coolftn(object): - def __init__(self, fname=os.path.join(dirpath, 'coolftn.txt')): + def __init__(self, fname=os.path.join(dirpath, "coolftn.txt")): cdf = np.loadtxt(fname) self.cool = cdf[:, 0] self.heat = cdf[:, 1] @@ -87,6 +82,7 @@ def get_heat(self, T1): return heat + # basic qunatities with renormalization @@ -95,7 +91,8 @@ def _ndensity(field, data): def _ram_pok_z(field, data): - return data["gas", "density"] * data["gas", "velocity_z"]**2 / kboltz + return data["gas", "density"] * data["gas", "velocity_z"] ** 2 / kboltz + # thermodynamics quantities @@ -132,13 +129,16 @@ def _dvelocity(field, data): def _dvelocity_mag(field, data): - return np.sqrt(data["gas", "velocity_x"]**2 + - data["gas", "dvelocity_y"]**2 + - data["gas", "velocity_z"]**2) + return np.sqrt( + data["gas", "velocity_x"] ** 2 + + data["gas", "dvelocity_y"] ** 2 + + data["gas", "velocity_z"] ** 2 + ) def _dkinetic_energy(field, data): - return 0.5 * data['gas', 'dvelocity_magnitude']**2 * data['gas', 'density'] + return 0.5 * data["gas", "dvelocity_magnitude"] ** 2 * data["gas", "density"] + # magnetic fields @@ -146,6 +146,7 @@ def _dkinetic_energy(field, data): def _mag_pok(field, data): return data["gas", "magnetic_pressure"] / kboltz + # chemistry @@ -158,9 +159,15 @@ def _C(field, data): def _H(field, data): - return (1. - (data["H2"] * 2 + 3 * data["H3+"] + - 2 * data["H2+"] + data["H+"] + data["HCO+"] + - data["CHx"] + data["OHx"])) + return 1.0 - ( + data["H2"] * 2 + + 3 * data["H3+"] + + 2 * data["H2+"] + + data["H+"] + + data["HCO+"] + + data["CHx"] + + data["OHx"] + ) def _O(field, data): @@ -188,27 +195,28 @@ def _tH2(field, data): def _electron(field, data): - return (data["He+"] + data["C+"] + data["HCO+"] + - data["H+"] + data["H3+"] + data["H2+"]) + return ( + data["He+"] + data["C+"] + data["HCO+"] + data["H+"] + data["H3+"] + data["H2+"] + ) def _temperature_chem(field, data): kb = 1.381e-16 - return yt.units.K * data["E"] / (1.5 * kb * - (1. - data["H2"] + data["e"] + GxHe)) + return yt.units.K * data["E"] / (1.5 * kb * (1.0 - data["H2"] + data["e"] + GxHe)) + # UV def _radiation_FUV(field, data): ds = data.ds - pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit)**2) + pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit) ** 2) return data["athena", "rad_energy_density1"] * pressure_unit def _radiation_EUV(field, data): ds = data.ds - pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit)**2) + pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit) ** 2) return data["athena", "rad_energy_density0"] * pressure_unit @@ -229,20 +237,23 @@ def _radiation_nHI(field, data): def _radiation_nHII(field, data): - return (data["gas", "density"] * - (1 - data["athena", "specific_scalar[0]"]) / (muH * mh)) + return ( + data["gas", "density"] * (1 - data["athena", "specific_scalar[0]"]) / (muH * mh) + ) def _radiation_ne(field, data): - return (data["gas", "density"] * - (1 - data["athena", "specific_scalar[0]"]) / (muH * mh)) + return ( + data["gas", "density"] * (1 - data["athena", "specific_scalar[0]"]) / (muH * mh) + ) + # basic def _pressure(field, data): ds = data.ds - pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit)**2) + pressure_unit = ds.mass_unit / (ds.length_unit * (ds.time_unit) ** 2) return data["press"] * pressure_unit @@ -250,102 +261,254 @@ def _density(field, data): return data["rho"] -def add_yt_fields(ds, chemistry=True, rad=False, - cooling=True, mhd=False, rotation=False, zdg=1.): - Gzdg = zdg - GxHe = 0.1 - GxC = 1.6e-4 * Gzdg - GxO = 3.2e-4 * Gzdg - GxSi = 1.7e-6 * Gzdg - ds.add_field(("gas", "nH"), function=_ndensity, - units='cm**(-3)', display_name=r'$n_{\rm H}$', sampling_type="cell") - ds.add_field(("gas", "H_nuclei_density"), function=_ndensity, - units='cm**(-3)', display_name=r'$n_{\rm H}$', sampling_type="cell") - ds.add_field(("gas", "ram_pok_z"), function=_ram_pok_z, - units='K*cm**(-3)', display_name=r'$P_{\rm turb}/k_{\rm B}$', - sampling_type="cell") - ds.add_field(("gas", "density"), function=_density, display_name="density", - sampling_type="cell") - ds.add_field(("gas", "pressure"), function=_pressure, units="g/(cm* s**2)", - display_name="pressure", sampling_type="cell") +def add_yt_fields( + ds, chemistry=True, rad=False, cooling=True, mhd=False, rotation=False, zdg=1.0 +): + # Gzdg = zdg + # GxHe = 0.1 + # GxC = 1.6e-4 * Gzdg + # GxO = 3.2e-4 * Gzdg + # GxSi = 1.7e-6 * Gzdg + ds.add_field( + ("gas", "nH"), + function=_ndensity, + units="cm**(-3)", + display_name=r"$n_{\rm H}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "H_nuclei_density"), + function=_ndensity, + units="cm**(-3)", + display_name=r"$n_{\rm H}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "ram_pok_z"), + function=_ram_pok_z, + units="K*cm**(-3)", + display_name=r"$P_{\rm turb}/k_{\rm B}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "density"), + function=_density, + display_name="density", + sampling_type="cell", + ) + ds.add_field( + ("gas", "pressure"), + function=_pressure, + units="g/(cm* s**2)", + display_name="pressure", + sampling_type="cell", + ) if chemistry: - ds.add_field(("gas", "Si"), function=_Si, - units="", display_name=r"${\rm Si}$", sampling_type="cell") - ds.add_field(("gas", "C"), function=_C, - units="", display_name=r"${\rm C}$", sampling_type="cell") - ds.add_field(("gas", "H"), function=_H, - units="", display_name=r"${\rm H}$", sampling_type="cell") - ds.add_field(("gas", "O"), function=_O, - units="", display_name=r"${\rm O}$", sampling_type="cell") - ds.add_field(("gas", "He"), function=_He, - units="", display_name=r"${\rm He}$", sampling_type="cell") - ds.add_field(("gas", "CO2Ctot"), function=_CO2Ctot, - units="", display_name=r"${\rm CO}/{\rm C}_{\rm tot}$", - sampling_type="cell") - ds.add_field(("gas", "C2Ctot"), function=_C2Ctot, - units="", display_name=r"${\rm C}/{\rm C}_{\rm tot}$", - sampling_type="cell") - ds.add_field(("gas", "C+2Ctot"), function=_Cp2Ctot, - units="", display_name=r"${\rm C^+}/{\rm C}_{\rm tot}$", - sampling_type="cell") - ds.add_field(("gas", "2H2"), function=_tH2, - units="", display_name=r"${\rm 2H_2}$", sampling_type="cell") - ds.add_field(("gas", "e"), function=_electron, - units="", display_name=r"${\rm e^{-}}$", sampling_type="cell") - ds.add_field(("gas", "temperature_chem"), function=_temperature_chem, - units="K", display_name=r"$T$", sampling_type="cell") + ds.add_field( + ("gas", "Si"), + function=_Si, + units="", + display_name=r"${\rm Si}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "C"), + function=_C, + units="", + display_name=r"${\rm C}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "H"), + function=_H, + units="", + display_name=r"${\rm H}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "O"), + function=_O, + units="", + display_name=r"${\rm O}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "He"), + function=_He, + units="", + display_name=r"${\rm He}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "CO2Ctot"), + function=_CO2Ctot, + units="", + display_name=r"${\rm CO}/{\rm C}_{\rm tot}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "C2Ctot"), + function=_C2Ctot, + units="", + display_name=r"${\rm C}/{\rm C}_{\rm tot}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "C+2Ctot"), + function=_Cp2Ctot, + units="", + display_name=r"${\rm C^+}/{\rm C}_{\rm tot}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "2H2"), + function=_tH2, + units="", + display_name=r"${\rm 2H_2}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "e"), + function=_electron, + units="", + display_name=r"${\rm e^{-}}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "temperature_chem"), + function=_temperature_chem, + units="K", + display_name=r"$T$", + sampling_type="cell", + ) if cooling: - ds.add_field(("gas", "pok"), function=_pok, - units='K*cm**(-3)', display_name=r'$P/k_{\rm B}$', - sampling_type="cell") - ds.add_field(("gas", "cs"), function=_cs, - units='km*s**(-1)', display_name=r'$c_s$', sampling_type="cell") - ds.add_field(("gas", "T1"), function=_T1, - units='K', display_name=r'$T_1$', sampling_type="cell") - ds.add_field(("gas", "mu_cgk"), function=_mu, - units='', display_name=r'$\mu_{\rm cgk}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "temperature"), function=_temperature, - units='K', display_name=r'$T$', force_override=True, - sampling_type="cell") + ds.add_field( + ("gas", "pok"), + function=_pok, + units="K*cm**(-3)", + display_name=r"$P/k_{\rm B}$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "cs"), + function=_cs, + units="km*s**(-1)", + display_name=r"$c_s$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "T1"), + function=_T1, + units="K", + display_name=r"$T_1$", + sampling_type="cell", + ) + ds.add_field( + ("gas", "mu_cgk"), + function=_mu, + units="", + display_name=r"$\mu_{\rm cgk}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "temperature"), + function=_temperature, + units="K", + display_name=r"$T$", + force_override=True, + sampling_type="cell", + ) if rotation: - ds.add_field(("gas", "dvelocity_y"), function=_dvelocity, - units='km/s', display_name=r'$\delta v_y$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "dvelocity_magnitude"), function=_dvelocity_mag, - units='km/s', display_name=r'$v$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "dkinetic_energy"), function=_dkinetic_energy, - units='erg/cm**3', display_name=r'$E_k$', force_override=True, - sampling_type="cell") + ds.add_field( + ("gas", "dvelocity_y"), + function=_dvelocity, + units="km/s", + display_name=r"$\delta v_y$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "dvelocity_magnitude"), + function=_dvelocity_mag, + units="km/s", + display_name=r"$v$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "dkinetic_energy"), + function=_dkinetic_energy, + units="erg/cm**3", + display_name=r"$E_k$", + force_override=True, + sampling_type="cell", + ) if mhd: - ds.add_field(("gas", "mag_pok"), function=_mag_pok, - units='K*cm**(-3)', display_name=r'$P_{\rm mag}/k_{\rm B}$', - sampling_type="cell") + ds.add_field( + ("gas", "mag_pok"), + function=_mag_pok, + units="K*cm**(-3)", + display_name=r"$P_{\rm mag}/k_{\rm B}$", + sampling_type="cell", + ) if rad: - ds.add_field(("gas", "nHI"), function=_radiation_nHI, - units='cm**(-3)', display_name=r'$n_{HI}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "xHI"), function=_radiation_xHI, - units='dimensionless', display_name=r'$x_{HI}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "nHII"), function=_radiation_nHII, - units='cm**(-3)', display_name=r'$n_{HII}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "xHII"), function=_radiation_xHII, - units='dimensionless', display_name=r'$x_{HII}$', - force_override=True, - sampling_type="cell") - ds.add_field(("gas", "ne"), function=_radiation_ne, - units='cm**(-3)', display_name=r'$n_{e}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "El_number_density"), function=_radiation_ne, - units='cm**(-3)', display_name=r'$n_{e}$', force_override=True, - sampling_type="cell") - ds.add_field(("gas", "xe"), function=_radiation_xe, - units='dimensionless', display_name=r'$x_{e}$', - force_override=True, - sampling_type="cell") + ds.add_field( + ("gas", "nHI"), + function=_radiation_nHI, + units="cm**(-3)", + display_name=r"$n_{HI}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "xHI"), + function=_radiation_xHI, + units="dimensionless", + display_name=r"$x_{HI}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "nHII"), + function=_radiation_nHII, + units="cm**(-3)", + display_name=r"$n_{HII}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "xHII"), + function=_radiation_xHII, + units="dimensionless", + display_name=r"$x_{HII}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "ne"), + function=_radiation_ne, + units="cm**(-3)", + display_name=r"$n_{e}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "El_number_density"), + function=_radiation_ne, + units="cm**(-3)", + display_name=r"$n_{e}$", + force_override=True, + sampling_type="cell", + ) + ds.add_field( + ("gas", "xe"), + function=_radiation_xe, + units="dimensionless", + display_name=r"$x_{e}$", + force_override=True, + sampling_type="cell", + ) # ds.add_field(("gas","Erad_EUV"), function=_radiation_EUV, \ # units='erg*cm**(-3)',display_name=r'$\mathcal{E}_{\rm EUV}$', @@ -357,58 +520,97 @@ def add_yt_fields(ds, chemistry=True, rad=False, # sampling_type="cell") -def set_aux(model='solar'): +def set_aux(model="solar"): aux = {} - aux['nH'] = dict(label=r'$n_H\;[{\rm cm}^{-3}]$', - unit='cm**(-3)', limits=(1.e-6, 1.e6), - cmap=plt.cm.Spectral_r, clim=(2.e-5, 2.e2), - cticks=(1.e-4, 1.e-2, 1, 1.e2), - n_bins=128, log=True) - aux['pok'] = dict(label=r'$P/k_B\;[{\rm K}\,{\rm cm}^{-3}]$', - unit='K*cm**(-3)', limits=(1.e-2, 1.e8), - cmap=plt.cm.gnuplot2, clim=(10, 5.e5), - n_bins=128, log=True) - aux['temperature'] = dict(label=r'$T\;[{\rm K}]$', - unit='K', limits=(1.e0, 1.e9), - cmap=plt.cm.RdYlBu_r, - clim=(10, 1.e8), - n_bins=128, log=True) - aux['surface_density'] = dict( - label=r'$\Sigma [{\rm M}_{\odot} {\rm pc}^{-2}]$', - cmap=plt.cm.pink_r, clim=(0.1, 100), log=True) - aux['dvelocity_magnitude'] = dict(label=r'$v [{\rm km/s}]$', - unit='km/s', limits=(0.1, 1.e4), - cmap=plt.cm.jet, clim=(1, 1000), - n_bins=128, log=True) - aux['velocity_z'] = dict(label=r'$v_z [{\rm km/s}]$', - unit='km/s', limits=(-1500, 1500), - cmap=plt.cm.RdBu_r, clim=(-200, 200), - cticks=(-100, 0, 100), - n_bins=256, log=False) - aux['magnetic_field_strength'] = dict(label=r'$B [\mu{\rm G}]$', - unit='gauss', - cmap=plt.cm.viridis, clim=(0.01, 10), - factor=1.e6, - n_bins=128, log=True) - aux['mag_pok'] = dict(label=r'$P_{\rm mag}/k_B [{\rm K}{\rm cm}^{-3}]$', - unit='K*cm**(-3)', limits=(1.e-2, 1.e8), - n_bins=128, log=True) - aux['ram_pok_z'] = dict( - label=r'$P_{\rm turb}/k_B [{\rm K}{\rm cm}^{-3}]$', - unit='K*cm**(-3)', limits=(1.e-2, 1.e8), - n_bins=128, log=True) - aux['plasma_beta'] = dict(label=r'$\beta$', limits=(1.e-4, 1.e16), - n_bins=256, log=True) - - if model == 'starburst': - aux['nH']['clim'] = (2.e-5, 2.e3) - aux['pok']['clim'] = (10, 5.e6) - aux['surface_density']['clim'] = (1, 1000) - aux['magnetic_field_strength']['clim'] = (0.1, 100) - - if model == 'multi_SN': - aux['nH']['clim'] = (2.e-5, 2.e2) - aux['pok']['clim'] = (50, 1.e5) + aux["nH"] = dict( + label=r"$n_H\;[{\rm cm}^{-3}]$", + unit="cm**(-3)", + limits=(1.0e-6, 1.0e6), + cmap=plt.cm.Spectral_r, + clim=(2.0e-5, 2.0e2), + cticks=(1.0e-4, 1.0e-2, 1, 1.0e2), + n_bins=128, + log=True, + ) + aux["pok"] = dict( + label=r"$P/k_B\;[{\rm K}\,{\rm cm}^{-3}]$", + unit="K*cm**(-3)", + limits=(1.0e-2, 1.0e8), + cmap=plt.cm.gnuplot2, + clim=(10, 5.0e5), + n_bins=128, + log=True, + ) + aux["temperature"] = dict( + label=r"$T\;[{\rm K}]$", + unit="K", + limits=(1.0e0, 1.0e9), + cmap=plt.cm.RdYlBu_r, + clim=(10, 1.0e8), + n_bins=128, + log=True, + ) + aux["surface_density"] = dict( + label=r"$\Sigma [{\rm M}_{\odot} {\rm pc}^{-2}]$", + cmap=plt.cm.pink_r, + clim=(0.1, 100), + log=True, + ) + aux["dvelocity_magnitude"] = dict( + label=r"$v [{\rm km/s}]$", + unit="km/s", + limits=(0.1, 1.0e4), + cmap=plt.cm.jet, + clim=(1, 1000), + n_bins=128, + log=True, + ) + aux["velocity_z"] = dict( + label=r"$v_z [{\rm km/s}]$", + unit="km/s", + limits=(-1500, 1500), + cmap=plt.cm.RdBu_r, + clim=(-200, 200), + cticks=(-100, 0, 100), + n_bins=256, + log=False, + ) + aux["magnetic_field_strength"] = dict( + label=r"$B [\mu{\rm G}]$", + unit="gauss", + cmap=plt.cm.viridis, + clim=(0.01, 10), + factor=1.0e6, + n_bins=128, + log=True, + ) + aux["mag_pok"] = dict( + label=r"$P_{\rm mag}/k_B [{\rm K}{\rm cm}^{-3}]$", + unit="K*cm**(-3)", + limits=(1.0e-2, 1.0e8), + n_bins=128, + log=True, + ) + aux["ram_pok_z"] = dict( + label=r"$P_{\rm turb}/k_B [{\rm K}{\rm cm}^{-3}]$", + unit="K*cm**(-3)", + limits=(1.0e-2, 1.0e8), + n_bins=128, + log=True, + ) + aux["plasma_beta"] = dict( + label=r"$\beta$", limits=(1.0e-4, 1.0e16), n_bins=256, log=True + ) + + if model == "starburst": + aux["nH"]["clim"] = (2.0e-5, 2.0e3) + aux["pok"]["clim"] = (10, 5.0e6) + aux["surface_density"]["clim"] = (1, 1000) + aux["magnetic_field_strength"]["clim"] = (0.1, 100) + + if model == "multi_SN": + aux["nH"]["clim"] = (2.0e-5, 2.0e2) + aux["pok"]["clim"] = (50, 1.0e5) return aux @@ -420,14 +622,18 @@ def check_aux(fields): print(aux[f]) -unit_base_pp = {"length_unit": (1.0, "pc"), - "time_unit": (1.0, "s*pc/km"), - "mass_unit": (2.38858753789e-24, "g/cm**3*pc**3")} -unit_base = {"length_unit": (1.0, "pc"), - "time_unit": (1.0, "s*pc/km"), - "mass_unit": (2.38858753789e-24, "g/cm**3*pc**3"), - "velocity_unit": (1.0, "km/s"), - "magnetic_unit": (5.4786746797e-07, "gauss")} +unit_base_pp = { + "length_unit": (1.0, "pc"), + "time_unit": (1.0, "s*pc/km"), + "mass_unit": (2.38858753789e-24, "g/cm**3*pc**3"), +} +unit_base = { + "length_unit": (1.0, "pc"), + "time_unit": (1.0, "s*pc/km"), + "mass_unit": (2.38858753789e-24, "g/cm**3*pc**3"), + "velocity_unit": (1.0, "km/s"), + "magnetic_unit": (5.4786746797e-07, "gauss"), +} def load_athena4p2_mhd(filename): @@ -458,25 +664,25 @@ def save_clumps(clumps_list, filename, fields_list=None): base_object = clumps_list[0].data.base_object if fields_list is None: fl = [ - ('gas', 'density'), - ('gas', 'magnetic_field_x'), - ('gas', 'magnetic_field_y'), - ('gas', 'magnetic_field_z'), - ('gas', 'pressure'), - ('gas', 'velocity_x'), - ('gas', 'velocity_y'), - ('gas', 'velocity_z'), - ('gravitational_potential'), - ('gas', 'phi'), - ('cell_mass'), - ('cell_volume') + ("gas", "density"), + ("gas", "magnetic_field_x"), + ("gas", "magnetic_field_y"), + ("gas", "magnetic_field_z"), + ("gas", "pressure"), + ("gas", "velocity_x"), + ("gas", "velocity_y"), + ("gas", "velocity_z"), + ("gravitational_potential"), + ("gas", "phi"), + ("cell_mass"), + ("cell_volume"), ] else: fl = copy.deepcopy(fields_list) conditional_output = [] # generate fields list for c in clumps_list: - contour_id = c.data.conditionals[0].split('contours')[1].split("'")[0] + contour_id = c.data.conditionals[0].split("contours")[1].split("'")[0] conditional_i = c.data.conditionals[0] fl.append(("index", "contours" + contour_id)) conditional_output.append(("contours" + contour_id, conditional_i)) @@ -495,62 +701,83 @@ def read_clumps(data_set, fn_cond=None): clumps = [] fo = open(fn_cond) for line in fo.readlines(): - ls = line.split(',') - cli = ls[0] - condi = ls[1].split('\n')[0] + ls = line.split(",") + # cli = ls[0] + condi = ls[1].split("\n")[0] print(condi) - cut = data_set.cut_region(data_set.all_data(), - [condi]) - c1 = Clump(cut, ('grid', 'phi')) + cut = data_set.cut_region(data_set.all_data(), [condi]) + c1 = Clump(cut, ("grid", "phi")) clumps.append(c1) return clumps def get_clump_properties(cl): - vol = cl['cell_volume'][0] - totvol = cl.quantities.total_quantity('cell_volume').in_units('pc**3') - radius = (totvol * 3. * math.pi / 4.)**(1. / 3.) - mass = cl.quantities.total_quantity('cell_mass').in_units('Msun') - nden = cl.quantities.weighted_average_quantity('density', 'cell_mass') - vx_c = cl.quantities.weighted_average_quantity( - "velocity_x", "cell_mass").in_units("km/s") - vy_c = cl.quantities.weighted_average_quantity( - "velocity_y", "cell_mass").in_units("km/s") - vz_c = cl.quantities.weighted_average_quantity( - "velocity_z", "cell_mass").in_units("km/s") + vol = cl["cell_volume"][0] + totvol = cl.quantities.total_quantity("cell_volume").in_units("pc**3") + radius = (totvol * 3.0 * math.pi / 4.0) ** (1.0 / 3.0) + mass = cl.quantities.total_quantity("cell_mass").in_units("Msun") + nden = cl.quantities.weighted_average_quantity("density", "cell_mass") + vx_c = cl.quantities.weighted_average_quantity("velocity_x", "cell_mass").in_units( + "km/s" + ) + vy_c = cl.quantities.weighted_average_quantity("velocity_y", "cell_mass").in_units( + "km/s" + ) + vz_c = cl.quantities.weighted_average_quantity("velocity_z", "cell_mass").in_units( + "km/s" + ) vc = [vx_c, vy_c, vz_c] - Ekin_cells = 0.5 * cl["cell_mass"] * ((cl["velocity_x"] - vx_c)**2 + - (cl["velocity_y"] - vy_c)**2 + - (cl["velocity_z"] - vz_c)**2) + Ekin_cells = ( + 0.5 + * cl["cell_mass"] + * ( + (cl["velocity_x"] - vx_c) ** 2 + + (cl["velocity_y"] - vy_c) ** 2 + + (cl["velocity_z"] - vz_c) ** 2 + ) + ) Ekin = Ekin_cells.sum() - Emag = cl.quantities.total_quantity('magnetic_energy') * vol.in_cgs() - Eth = cl.quantities.total_quantity('pressure') * vol.in_cgs() * 1.5 - Egrav = -(cl["gravitational_potential"].max() - - cl["gravitational_potential"].min()) * (u.km / u.s)**2 * mass - B = yt.YTArray([cl.quantities.weighted_average_quantity('magnetic_field_x', - 'cell_mass'), - cl.quantities.weighted_average_quantity( - 'magnetic_field_y', 'cell_mass'), - cl.quantities.weighted_average_quantity('magnetic_field_z', - 'cell_mass')]) + Emag = cl.quantities.total_quantity("magnetic_energy") * vol.in_cgs() + Eth = cl.quantities.total_quantity("pressure") * vol.in_cgs() * 1.5 + Egrav = ( + -(cl["gravitational_potential"].max() - cl["gravitational_potential"].min()) + * (u.km / u.s) ** 2 + * mass + ) + B = yt.YTArray( + [ + cl.quantities.weighted_average_quantity("magnetic_field_x", "cell_mass"), + cl.quantities.weighted_average_quantity("magnetic_field_y", "cell_mass"), + cl.quantities.weighted_average_quantity("magnetic_field_z", "cell_mass"), + ] + ) Eratio = -((Ekin + Emag + Eth) / Egrav) - print('==========================================') - print('Volume:', totvol) - print('Effective radius:', radius) - print('Mass:', mass) - print('center of mass velocity', vc) - print('Mass-weighted mean Density:', nden / u.mh) - print('Mass-weighted mean B (micro G):', B.in_units('uG')) - print('Ekin:', Ekin.in_units('erg')) - print('Emag:', Emag.in_units('erg')) - print('Eth:', Eth.in_units('erg')) - print('Egrav:', Egrav.in_units('erg')) - print('Eratio = (Ekin+Emag+Eth)/(-Egrav):', Eratio.in_cgs()) - print('==========================================') - return {"volume": totvol, "radius": radius, "mtot": mass, "vc": vc, - "density_mass_weighted": nden, "B_mass_weigthed": B, - "Ekin": Ekin, "Emag": Emag, "Eth": Eth, - "Egrav": Egrav, "Eratio": Eratio} + print("==========================================") + print("Volume:", totvol) + print("Effective radius:", radius) + print("Mass:", mass) + print("center of mass velocity", vc) + print("Mass-weighted mean Density:", nden / u.mh) + print("Mass-weighted mean B (micro G):", B.in_units("uG")) + print("Ekin:", Ekin.in_units("erg")) + print("Emag:", Emag.in_units("erg")) + print("Eth:", Eth.in_units("erg")) + print("Egrav:", Egrav.in_units("erg")) + print("Eratio = (Ekin+Emag+Eth)/(-Egrav):", Eratio.in_cgs()) + print("==========================================") + return { + "volume": totvol, + "radius": radius, + "mtot": mass, + "vc": vc, + "density_mass_weighted": nden, + "B_mass_weigthed": B, + "Ekin": Ekin, + "Emag": Emag, + "Eth": Eth, + "Egrav": Egrav, + "Eratio": Eratio, + } def get_covering_grid(ds, left_edge=None, dims=None): @@ -559,8 +786,7 @@ def get_covering_grid(ds, left_edge=None, dims=None): left_edge = ds.domain_left_edge if dims is None: dims = ds.domain_dimensions - return ds.covering_grid(level=0, left_edge=left_edge, - dims=dims) + return ds.covering_grid(level=0, left_edge=left_edge, dims=dims) def get_extent(grid, axis, unit): diff --git a/miscellaneous/read_vtk.py b/miscellaneous/read_vtk.py index 222e9e6..40e4d18 100644 --- a/miscellaneous/read_vtk.py +++ b/miscellaneous/read_vtk.py @@ -6,14 +6,14 @@ import os import os.path as osp -import glob, struct +import glob +import struct import numpy as np import xarray as xr -import astropy.constants as ac -import astropy.units as au from units import Units + def read_vtk(filename, id0_only=False): """Convenience wrapper function to read Athena vtk output file using AthenaDataSet class. @@ -33,8 +33,8 @@ def read_vtk(filename, id0_only=False): return AthenaDataSet(filename, id0_only=id0_only) -class AthenaDataSet(object): +class AthenaDataSet(object): def __init__(self, filename, id0_only=False, units=Units(), dfi=None): """Class to read athena vtk file. @@ -52,10 +52,11 @@ def __init__(self, filename, id0_only=False, units=Units(), dfi=None): """ if not osp.exists(filename): - raise IOError(('File does not exist: {0:s}'.format(filename))) + raise IOError(("File does not exist: {0:s}".format(filename))) - dirname, problem_id, num, suffix, ext, mpi_mode, nonzero_id = \ - _parse_filename(filename) + dirname, problem_id, num, suffix, ext, mpi_mode, nonzero_id = _parse_filename( + filename + ) if id0_only: mpi_mode = False @@ -77,16 +78,22 @@ def __init__(self, filename, id0_only=False, units=Units(), dfi=None): # Find all vtk file names and add to flist if mpi_mode: if self.suffix is None: - fname_pattern = osp.join(dirname, 'id*/{0:s}-id*.{1:s}.{2:s}'.\ - format(problem_id, num, ext)) + fname_pattern = osp.join( + dirname, "id*/{0:s}-id*.{1:s}.{2:s}".format(problem_id, num, ext) + ) else: - fname_pattern = osp.join(dirname, 'id*/{0:s}-id*.{1:s}.{2:s}.{3:s}'.\ - format(problem_id, num, suffix, ext)) + fname_pattern = osp.join( + dirname, + "id*/{0:s}-id*.{1:s}.{2:s}.{3:s}".format( + problem_id, num, suffix, ext + ), + ) fnames = glob.glob(fname_pattern) self.fnames += fnames if nonzero_id: from collections import OrderedDict + self.fnames = list(OrderedDict.fromkeys(self.fnames)) self.grid = self._set_grid() @@ -94,18 +101,17 @@ def __init__(self, filename, id0_only=False, units=Units(), dfi=None): self.set_region() # Need separte field_map for different grids - if self.domain['all_grid_equal']: + if self.domain["all_grid_equal"]: self._field_map = _set_field_map(self.grid[0]) for g in self.grid: - g['field_map'] = self._field_map + g["field_map"] = self._field_map else: for g in self.grid: - g['field_map'] = _set_field_map(g) - self._field_map = self.grid[0]['field_map'] + g["field_map"] = _set_field_map(g) + self._field_map = self.grid[0]["field_map"] self.field_list = list(self._field_map.keys()) - def get_cc_pos(self): """Compute cell center positions @@ -116,9 +122,10 @@ def get_cc_pos(self): """ xc = dict() - for axis, le, re, dx in zip(('x', 'y', 'z'), \ - self.region['gle'], self.region['gre'], self.domain['dx']): - xc[axis] = np.arange(le + 0.5*dx, re + 0.5*dx, dx) + for axis, le, re, dx in zip( + ("x", "y", "z"), self.region["gle"], self.region["gre"], self.domain["dx"] + ): + xc[axis] = np.arange(le + 0.5 * dx, re + 0.5 * dx, dx) return xc @@ -131,44 +138,50 @@ def get_cc_ijk(s, x1, x2, x3): Particle position """ domain = s.domain - le1,le2,le3 = domain['le'] - dx1,dx2,dx3 = domain['dx'] + le1, le2, le3 = domain["le"] + dx1, dx2, dx3 = domain["dx"] - return (np.floor((x1 - le1)/dx1).astype(int), - np.floor((x2 - le2)/dx2).astype(int), - np.floor((x3 - le3)/dx3).astype(int)) + return ( + np.floor((x1 - le1) / dx1).astype(int), + np.floor((x2 - le2) / dx2).astype(int), + np.floor((x3 - le3) / dx3).astype(int), + ) def set_region(self, le=None, re=None): - """Set region and find overlapping grids. - """ + """Set region and find overlapping grids.""" if le is None: - le = self.domain['le'] + le = self.domain["le"] if re is None: - re = self.domain['re'] + re = self.domain["re"] le = np.array(le) re = np.array(re) if (re < le).any(): - raise ValueError('Check left/right edge.') + raise ValueError("Check left/right edge.") # Find all overlapping grids and their edges gle_all = [] # grid left edge gre_all = [] # grid right edge - gidx = [] # grid indices that belongs to this region - #print(self.grid,len(self.grid)) + gidx = [] # grid indices that belongs to this region + # print(self.grid,len(self.grid)) for i, g in enumerate(self.grid): - if (g['re'] >= le).all() and (g['le'] <= re).all(): + if (g["re"] >= le).all() and (g["le"] <= re).all(): gidx.append(i) - gle_all.append(g['le']) - gre_all.append(g['re']) + gle_all.append(g["le"]) + gre_all.append(g["re"]) gidx = np.array(gidx) if len(gidx) == 0: - raise ValueError('Check left/right edges:', le, re, \ - ' Domain left/right edges are ', \ - self.domain['le'], self.domain['re']) + raise ValueError( + "Check left/right edges:", + le, + re, + " Domain left/right edges are ", + self.domain["le"], + self.domain["re"], + ) gle_all = np.array(gle_all) gre_all = np.array(gre_all) @@ -185,10 +198,13 @@ def set_region(self, le=None, re=None): NGrid = np.array([len(gleu_) for gleu_ in gleu]) # Number of cells per grid - Nxg = np.concatenate([(greu_ - gleu_)/dx for greu_, gleu_, dx in \ - zip(greu, gleu, self.domain['dx'])]) - Nxg = np.array(np.array_split(Nxg, NGrid.cumsum()[:-1]), - dtype=object).tolist() + Nxg = np.concatenate( + [ + (greu_ - gleu_) / dx + for greu_, gleu_, dx in zip(greu, gleu, self.domain["dx"]) + ] + ) + Nxg = np.array(np.array_split(Nxg, NGrid.cumsum()[:-1]), dtype=object).tolist() # Since floating point arithmetic may result in incorrect results, need # to round to the nearest integer @@ -200,17 +216,26 @@ def set_region(self, le=None, re=None): for i, Nxg_ in enumerate(Nxg): Nxr[i] = np.sum(Nxg_) - #print(gidx,NGrid) - assert len(gidx) == NGrid.prod(),\ - print('Unexpected error: Number of grids {0:d} != '.format(len(gidx)) + - 'number of unique edges {0:d}.'.format(NGrid.prod())) - - self.region = dict(le=le, re=re, gidx=gidx, - gleu=gleu, greu=greu,\ - gle=gle, gre=gre, - NGrid=NGrid, Nxg=Nxg, Nxr=Nxr) - - def get_slice(self, axis, field='density', pos='c', method='nearest'): + # print(gidx,NGrid) + assert len(gidx) == NGrid.prod(), print( + "Unexpected error: Number of grids {0:d} != ".format(len(gidx)) + + "number of unique edges {0:d}.".format(NGrid.prod()) + ) + + self.region = dict( + le=le, + re=re, + gidx=gidx, + gleu=gleu, + greu=greu, + gle=gle, + gre=gre, + NGrid=NGrid, + Nxg=Nxg, + Nxr=Nxr, + ) + + def get_slice(self, axis, field="density", pos="c", method="nearest"): """Read slice of fields. Parameters @@ -234,28 +259,28 @@ def get_slice(self, axis, field='density', pos='c', method='nearest'): axis_idx = dict(x=0, y=1, z=2) if pos is None: - pos = 'c' + pos = "c" field = np.atleast_1d(field) axis = np.atleast_1d(axis) for ax in axis: - le = np.copy(self.domain['le']) - re = np.copy(self.domain['re']) - if pos in ['c', 'center']: - pos = self.domain['center'][axis_idx[ax]] + le = np.copy(self.domain["le"]) + re = np.copy(self.domain["re"]) + if pos in ["c", "center"]: + pos = self.domain["center"][axis_idx[ax]] # Let's make sure le < re always and truncation error does not cause # problem, although performance can be slowed down a bit. - le[axis_idx[ax]] = pos - 0.5*self.domain['dx'][axis_idx[ax]] - re[axis_idx[ax]] = pos + 0.5*self.domain['dx'][axis_idx[ax]] + le[axis_idx[ax]] = pos - 0.5 * self.domain["dx"][axis_idx[ax]] + re[axis_idx[ax]] = pos + 0.5 * self.domain["dx"][axis_idx[ax]] dat = self.get_field(field, le, re, as_xarray=True) - slc = dat.sel(method='nearest', **{ax:pos}) + slc = dat.sel(method="nearest", **{ax: pos}) return slc - def get_field(self, field='density', le=None, re=None, as_xarray=True): + def get_field(self, field="density", le=None, re=None, as_xarray=True): """Read 3d fields data. Parameters @@ -293,7 +318,7 @@ def get_field(self, field='density', le=None, re=None, as_xarray=True): if not dflist.issubset(set(self.dfi.keys())): tmp = [] for f in dflist: - if not f in self.dfi.keys(): + if f not in self.dfi.keys(): tmp.append(f) raise KeyError("Unrecognized field name(s):", tmp) @@ -304,16 +329,16 @@ def get_field(self, field='density', le=None, re=None, as_xarray=True): # Fields that need to be read to calculate derived field flist_dep = set() for f in dflist: - flist_dep = flist_dep | set(self.dfi[f]['field_dep']) + flist_dep = flist_dep | set(self.dfi[f]["field_dep"]) # Fields names to be dropped later fdrop_list = flist_dep - flist # Need to adjust names for vector fields for f in list(fdrop_list): - if self._field_map[f]['nvar'] > 1: - for i in range(self._field_map[f]['nvar']): - fdrop_list.add(f+str(i+1)) + if self._field_map[f]["nvar"] > 1: + for i in range(self._field_map[f]["nvar"]): + fdrop_list.add(f + str(i + 1)) fdrop_list.remove(f) field = list(flist_dep | flist) @@ -321,20 +346,19 @@ def get_field(self, field='density', le=None, re=None, as_xarray=True): # Calculate derived fields for f in dflist: - dat[f] = self.dfi[f]['func'](dat, self.u) + dat[f] = self.dfi[f]["func"](dat, self.u) # Drop fields that are not requested if as_xarray: dat = dat.drop(list(fdrop_list)) - dat.attrs['dfi'] = self.dfi + dat.attrs["dfi"] = self.dfi else: for f in fdrop_list: del dat[f] return dat.squeeze() - def _get_field(self, field='density', le=None, re=None, as_xarray=True): - + def _get_field(self, field="density", le=None, re=None, as_xarray=True): field = np.atleast_1d(field) # Create region @@ -345,19 +369,23 @@ def _get_field(self, field='density', le=None, re=None, as_xarray=True): if as_xarray: # Cell center positions x = dict() - for axis, le, re, dx in zip(('x', 'y', 'z'), \ - self.region['gle'], self.region['gre'], self.domain['dx']): + for axis, le, re, dx in zip( + ("x", "y", "z"), + self.region["gle"], + self.region["gre"], + self.domain["dx"], + ): # May not result in correct number of elements due to truncation error # x[axis] = np.arange(le + 0.5*dx, re + 0.5*dx, dx) - x[axis] = np.arange(le + 0.5*dx, re + 0.25*dx, dx) + x[axis] = np.arange(le + 0.5 * dx, re + 0.25 * dx, dx) dat = dict() for k, v in arr.items(): - if len(v.shape) > self.domain['ndim']: + if len(v.shape) > self.domain["ndim"]: for i in range(v.shape[-1]): - dat[k + str(i+1)] = (('z','y','x'), v[..., i]) + dat[k + str(i + 1)] = (("z", "y", "x"), v[..., i]) else: - dat[k] = (('z','y','x'), v) + dat[k] = (("z", "y", "x"), v) attrs = dict() for k, v in self.domain.items(): @@ -371,19 +399,18 @@ def _get_field(self, field='density', le=None, re=None, as_xarray=True): return arr def _get_array(self, field): - arr = dict() for f in field: arr[f] = self._set_array(f) # Read from individual grids and copy to data - le = self.region['gle'] - dx = self.domain['dx'] - for i in self.region['gidx']: + le = self.region["gle"] + dx = self.domain["dx"] + for i in self.region["gidx"]: g = self.grid[i] - il = (np.rint((g['le'] - le)/dx)).astype(int) - iu = il + g['Nx'] - slc = tuple([slice(l, u) for l, u in zip(il[::-1], iu[::-1])]) + il = (np.rint((g["le"] - le) / dx)).astype(int) + iu = il + g["Nx"] + slc = tuple([slice(l_, u_) for l_, u_ in zip(il[::-1], iu[::-1])]) for f in field: arr[f][slc] = self._read_array(g, f) @@ -391,41 +418,38 @@ def _get_array(self, field): return arr def _read_array(self, grid, field): - - if field in grid['data']: - return grid['data'][field] + if field in grid["data"]: + return grid["data"][field] elif field in self.field_list: - fm = grid['field_map'][field] - if 'tarinfo' in grid: - fp = self.tarfile.extractfile(grid['tarinfo']) + fm = grid["field_map"][field] + if "tarinfo" in grid: + fp = self.tarfile.extractfile(grid["tarinfo"]) else: - fp = open(grid['filename'], 'rb') - fp.seek(fm['offset']) - fp.readline() # skip header - if fm['read_table']: + fp = open(grid["filename"], "rb") + fp.seek(fm["offset"]) + fp.readline() # skip header + if fm["read_table"]: fp.readline() - grid['data'][field] = np.asarray( - struct.unpack('>' + fm['ndata']*fm['dtype'], - fp.read(fm['dsize']))) + grid["data"][field] = np.asarray( + struct.unpack(">" + fm["ndata"] * fm["dtype"], fp.read(fm["dsize"])) + ) fp.close() - if fm['nvar'] == 1: - shape = np.flipud(grid['Nx']) + if fm["nvar"] == 1: + shape = np.flipud(grid["Nx"]) else: - shape = (*np.flipud(grid['Nx']), fm['nvar']) - - grid['data'][field].shape = shape + shape = (*np.flipud(grid["Nx"]), fm["nvar"]) - return grid['data'][field] + grid["data"][field].shape = shape + return grid["data"][field] def _set_array(self, field): - - dtype = self._field_map[field]['dtype'] - nvar = self._field_map[field]['nvar'] - Nxr = self.region['Nxr'] - if 'face_centered_B' in field: - Nxr[int(field[-1])-1] += 1 + dtype = self._field_map[field]["dtype"] + nvar = self._field_map[field]["nvar"] + Nxr = self.region["Nxr"] + if "face_centered_B" in field: + Nxr[int(field[-1]) - 1] += 1 if nvar == 1: shape = np.flipud(Nxr) else: @@ -433,42 +457,40 @@ def _set_array(self, field): return np.empty(shape, dtype=dtype) - def _set_domain(self): - domain = dict() grid = self.grid ngrid = len(grid) # Grid left/right edges - gle = np.empty((ngrid, 3), dtype='float32') - gre = np.empty((ngrid, 3), dtype='float32') - dx = np.empty((ngrid, 3), dtype='float32') - Nx = np.ones_like(dx, dtype='int') + gle = np.empty((ngrid, 3), dtype="float32") + gre = np.empty((ngrid, 3), dtype="float32") + dx = np.empty((ngrid, 3), dtype="float32") + Nx = np.ones_like(dx, dtype="int") for i, g in enumerate(grid): - gle[i, :] = g['le'] - gre[i, :] = g['re'] - Nx[i, :] = g['Nx'] - dx[i, :] = g['dx'] + gle[i, :] = g["le"] + gre[i, :] = g["re"] + Nx[i, :] = g["Nx"] + dx[i, :] = g["dx"] # Check if all grids have the equal size if (Nx[0] == Nx).all(): - domain['all_grid_equal'] = True + domain["all_grid_equal"] = True else: - domain['all_grid_equal'] = False + domain["all_grid_equal"] = False # Set domain le = gle.min(axis=0) re = gre.max(axis=0) - domain['ngrid'] = ngrid - domain['le'] = le - domain['re'] = re - domain['dx'] = dx[0, :] - domain['Lx'] = re - le - domain['center'] = 0.5*(le + re) - domain['Nx'] = np.round(domain['Lx']/domain['dx']).astype('int') - domain['ndim'] = 3 # always 3d + domain["ngrid"] = ngrid + domain["le"] = le + domain["re"] = re + domain["dx"] = dx[0, :] + domain["Lx"] = re - le + domain["center"] = 0.5 * (le + re) + domain["Nx"] = np.round(domain["Lx"] / domain["dx"]).astype("int") + domain["ndim"] = 3 # always 3d # file = open(self.fnames[0], 'rb') # tmpgrid = dict() @@ -479,7 +501,7 @@ def _set_domain(self): # file.close() # domain['time'] = tmpgrid['time'] - domain['time'] = grid[0]['time'] + domain["time"] = grid[0]["time"] return domain @@ -487,31 +509,30 @@ def _set_grid(self): grid = [] # Record filename and data_offset for i, fname in enumerate(self.fnames): - file = open(fname, 'rb') + file = open(fname, "rb") g = dict() - g['data'] = dict() - g['filename'] = fname - g['read_field'] = None - g['read_type'] = None + g["data"] = dict() + g["filename"] = fname + g["read_field"] = None + g["read_type"] = None - while g['read_field'] is None: - g['data_offset'] = file.tell() + while g["read_field"] is None: + g["data_offset"] = file.tell() line = file.readline() _vtk_parse_line(line, g) file.close() - g['Nx'] -= 1 - g['Nx'][g['Nx'] == 0] = 1 - g['dx'][g['Nx'] == 1] = 1.0 + g["Nx"] -= 1 + g["Nx"][g["Nx"] == 0] = 1 + g["dx"][g["Nx"] == 1] = 1.0 # Right edge - g['re'] = g['le'] + g['Nx']*g['dx'] + g["re"] = g["le"] + g["Nx"] * g["dx"] grid.append(g) return grid - def _parse_filename(filename): """Break up a filename into its component to check the extension and extract the output number. @@ -542,24 +563,23 @@ def _parse_filename(filename): nonzero_id = False # Check if dirname ends with id0 dirname_last = dirname.split(sep)[-1] - if dirname_last.startswith('id') and \ - dirname_last[2:].isdigit(): + if dirname_last.startswith("id") and dirname_last[2:].isdigit(): dirname = sep.join(dirname.split(sep)[:-1]) mpi_mode = True else: mpi_mode = False base = os.path.basename(filename) - base_split = base.split('.') + base_split = base.split(".") if len(base_split) == 3: - problem_id = '.'.join(base_split[:-2]) + problem_id = ".".join(base_split[:-2]) num = base_split[-2] suffix = None ext = base_split[-1] else: try: inum = -3 - test = int(base_split[inum]) + # test = int(base_split[inum]) # If dirname is idXX where XX>0, (2d vtk slices) # need to remove idXX string from the problem_id suffix = base_split[-2] @@ -568,11 +588,11 @@ def _parse_filename(filename): suffix = None if mpi_mode and int(dirname_last[2:]) != 0: - problem_id = '.'.join(base_split[:inum]) - problem_id = problem_id.replace('-' + dirname_last,'') + problem_id = ".".join(base_split[:inum]) + problem_id = problem_id.replace("-" + dirname_last, "") nonzero_id = True else: - problem_id = '.'.join(base_split[:inum]) + problem_id = ".".join(base_split[:inum]) num = base_split[inum] ext = base_split[-1] @@ -580,54 +600,52 @@ def _parse_filename(filename): return dirname, problem_id, num, suffix, ext, mpi_mode, nonzero_id - def _set_field_map(grid): - - fp = open(grid['filename'], 'rb') + fp = open(grid["filename"], "rb") fp.seek(0, 2) eof = fp.tell() - offset = grid['data_offset'] + offset = grid["data_offset"] fp.seek(offset) field_map = dict() - if 'Nx' in grid: - Nx = grid['Nx'] + if "Nx" in grid: + Nx = grid["Nx"] while offset < eof: line = fp.readline() sp = line.strip().split() - field = sp[1].decode('utf-8') + field = sp[1].decode("utf-8") field_map[field] = dict() - field_map[field]['read_table'] = False - field_map[field]['title'] = line + field_map[field]["read_table"] = False + field_map[field]["title"] = line if b"SCALARS" in line: tmp = fp.readline() - field_map[field]['read_table'] = True - field_map[field]['nvar'] = 1 + field_map[field]["read_table"] = True + field_map[field]["nvar"] = 1 elif b"VECTORS" in line: - field_map[field]['nvar'] = 3 + field_map[field]["nvar"] = 3 else: - raise TypeError(sp[0] + ' is unknown type.') - - field_map[field]['offset'] = offset - field_map[field]['ndata'] = field_map[field]['nvar']*grid['ncells'] - if field == 'face_centered_B1': - field_map[field]['ndata'] = (Nx[0]+1)*Nx[1]*Nx[2] - elif field == 'face_centered_B2': - field_map[field]['ndata'] = Nx[0]*(Nx[1]+1)*Nx[2] - elif field == 'face_centered_B3': - field_map[field]['ndata'] = Nx[0]*Nx[1]*(Nx[2]+1) - - if sp[2] == b'int': - dtype = 'i' - elif sp[2] == b'float': - dtype = 'f' - elif sp[2] == b'double': - dtype = 'd' - - field_map[field]['dtype'] = dtype - field_map[field]['dsize'] = field_map[field]['ndata']*struct.calcsize(dtype) - fp.seek(field_map[field]['dsize'], 1) + raise TypeError(sp[0] + " is unknown type.") + + field_map[field]["offset"] = offset + field_map[field]["ndata"] = field_map[field]["nvar"] * grid["ncells"] + if field == "face_centered_B1": + field_map[field]["ndata"] = (Nx[0] + 1) * Nx[1] * Nx[2] + elif field == "face_centered_B2": + field_map[field]["ndata"] = Nx[0] * (Nx[1] + 1) * Nx[2] + elif field == "face_centered_B3": + field_map[field]["ndata"] = Nx[0] * Nx[1] * (Nx[2] + 1) + + if sp[2] == b"int": + dtype = "i" + elif sp[2] == b"float": + dtype = "f" + elif sp[2] == b"double": + dtype = "d" + + field_map[field]["dtype"] = dtype + field_map[field]["dsize"] = field_map[field]["ndata"] * struct.calcsize(dtype) + fp.seek(field_map[field]["dsize"], 1) offset = fp.tell() tmp = fp.readline() if len(tmp) > 1: @@ -642,29 +660,29 @@ def _vtk_parse_line(line, grid): sp = line.strip().split() if b"vtk" in sp: - grid['vtk_version'] = sp[-1] + grid["vtk_version"] = sp[-1] elif b"time=" in sp: time_index = sp.index(b"time=") - grid['time'] = float(sp[time_index + 1].rstrip(b',')) - if b'level' in sp: - grid['level'] = int(sp[time_index + 3].rstrip(b',')) - if b'domain' in sp: - grid['domain'] = int(sp[time_index + 5].rstrip(b',')) + grid["time"] = float(sp[time_index + 1].rstrip(b",")) + if b"level" in sp: + grid["level"] = int(sp[time_index + 3].rstrip(b",")) + if b"domain" in sp: + grid["domain"] = int(sp[time_index + 5].rstrip(b",")) if sp[0] == b"PRIMITIVE": - grid['prim_var_type'] = True + grid["prim_var_type"] = True elif b"DIMENSIONS" in sp: - grid['Nx'] = np.array(sp[-3:]).astype('int') - elif b"ORIGIN" in sp: # left_edge - grid['le'] = np.array(sp[-3:]).astype('float64') + grid["Nx"] = np.array(sp[-3:]).astype("int") + elif b"ORIGIN" in sp: # left_edge + grid["le"] = np.array(sp[-3:]).astype("float64") elif b"SPACING" in sp: - grid['dx'] = np.array(sp[-3:]).astype('float64') + grid["dx"] = np.array(sp[-3:]).astype("float64") elif b"CELL_DATA" in sp: - grid['ncells'] = int(sp[-1]) + grid["ncells"] = int(sp[-1]) elif b"SCALARS" in sp: - grid['read_field'] = sp[1] - grid['read_type'] = 'scalar' + grid["read_field"] = sp[1] + grid["read_type"] = "scalar" elif b"VECTORS" in sp: - grid['read_field'] = sp[1] - grid['read_type'] = 'vector' + grid["read_field"] = sp[1] + grid["read_type"] = "vector" elif b"POINTS" in sp: - grid['npoint'] = eval(sp[1]) + grid["npoint"] = eval(sp[1]) diff --git a/prepare_tutorial.py b/prepare_tutorial.py index e0353c1..76f8b6c 100644 --- a/prepare_tutorial.py +++ b/prepare_tutorial.py @@ -1,6 +1,6 @@ -import sys -sys.path.insert(0, "./") # add path for import -import tigress_read # scripts for reading data +# import sys +# sys.path.insert(0, "./") # add path for import +import tigress_read # scripts for reading data # Master directory where the data is stored dir_master = "./data/" @@ -14,9 +14,9 @@ ## R8_4pc MHD while True: download_ok = input("Download 4pc MHD data [750MB]? (y/n):") - if download_ok.lower() in ['y','n']: + if download_ok.lower() in ["y", "n"]: break -if download_ok.lower() == 'y': +if download_ok.lower() == "y": # read the model information model_4pc = tigress_read.Model("R8_4pc", dir_master=dir_master) # download data @@ -27,9 +27,9 @@ ## R8_4pc MHD_PI while True: download_ok = input("Download 4pc MHD_PI data [3.9GB]? (y/n):") - if download_ok.lower() in ['y','n']: + if download_ok.lower() in ["y", "n"]: break -if download_ok.lower() == 'y': +if download_ok.lower() == "y": # read the model information model_4pc = tigress_read.Model("R8_4pc", dir_master=dir_master) # download data @@ -40,22 +40,22 @@ while True: download_ok = input("Download 2pc MHD/CO data [6GB]? (y/n):") - if download_ok.lower() in ['y','n']: + if download_ok.lower() in ["y", "n"]: break -if download_ok.lower() == 'y': - model_id_2pc = "R8_2pc" #name of the simulation model +if download_ok.lower() == "y": + model_id_2pc = "R8_2pc" # name of the simulation model model_2pc = tigress_read.Model(model_id_2pc, dir_master=dir_master) - #reading the model information + # reading the model information model_2pc.download(num, dataset="history") model_2pc.download(num, dataset="input") model_2pc.download(num, dataset="MHD") - model_2pc.download(num, dataset="CO_lines", iline=1) #CO(J=1-0) - model_2pc.download(num, dataset="CO_lines", iline=2) #CO(J=2-1) + model_2pc.download(num, dataset="CO_lines", iline=1) # CO(J=1-0) + model_2pc.download(num, dataset="CO_lines", iline=2) # CO(J=2-1) while True: download_ok = input("Download 2pc full chemistry data [15GB]? (y/n):") - if download_ok.lower() in ['y','n']: + if download_ok.lower() in ["y", "n"]: break -if download_ok.lower() == 'y': - model_2pc.download(num, dataset="chem") #large file 15GB +if download_ok.lower() == "y": + model_2pc.download(num, dataset="chem") # large file 15GB