diff --git a/csalt/csalt_disk.py b/csalt/csalt_disk.py index fc2425a..0c90898 100644 --- a/csalt/csalt_disk.py +++ b/csalt/csalt_disk.py @@ -3,7 +3,7 @@ from scipy.spatial import Delaunay import scipy.constants as sc import warnings - +from scipy.interpolate import LinearNDInterpolator class csalt_disk: @@ -85,6 +85,8 @@ def deproj_coords(self, dx, dy, incl, PA, z_func=None, else: y_dep = np.maximum.accumulate(y_dep, axis=0) + + """ # interpolation points_2D = np.column_stack((x_dep.flatten(), y_dep.flatten())) xsys = np.column_stack((xf.flatten(), yf.flatten() * np.cos(inc))) @@ -97,12 +99,23 @@ def deproj_coords(self, dx, dy, incl, PA, z_func=None, t_obs = np.arctan2(y_obs, az_sgn * x_obs) z_obs = z_func(r_obs) - # backfill extrapolations to a flat disk - r = np.where(np.isfinite(r_obs), r_obs, rf) - t = np.where(np.isfinite(t_obs), t_obs, tf) - z = np.where(np.isfinite(z_obs), z_obs, zf) + """ + #replacing interpolation with SciPy LinearNDInterpolator because the old interpolation method was leaving random NaNs in the middle of the grid + + interp_x = LinearNDInterpolator(list(zip(x_dep.flatten(), y_dep.flatten())), x_disk.flatten()) + interp_y = LinearNDInterpolator(list(zip(x_dep.flatten(), y_dep.flatten())), y_disk.flatten()) + + x_obs = np.reshape(interp_x(list(zip(xf.flatten(), yf.flatten()* np.cos(inc)))), self.x_sky.shape) + y_obs = np.reshape(interp_y(list(zip(xf.flatten(), yf.flatten()* np.cos(inc)))), self.y_sky.shape) + + r_obs = np.hypot(x_obs, y_obs) + t_obs = np.arctan2(y_obs, az_sgn * x_obs) + z_obs = z_func(r_obs) + + return r_obs, t_obs, z_obs + - return r, t, z + def set_disk_coords(self, dx, dy, incl, PA, dist, z_func=None, @@ -224,12 +237,12 @@ def set_vtheta_profile(self, function, side='both', vmin=0.0, vmax=None): @property def vtheta_f_proj(self): vt = self.vtheta_f * np.cos(self.t_disk_f) - return vt * np.sin(abs(np.radians(self.incl_f))) + return np.nan_to_num(vt * np.sin(abs(np.radians(self.incl_f)))) @property def vtheta_b_proj(self): vt = self.vtheta_b * np.cos(self.t_disk_b) - return vt * np.sin(abs(np.radians(self.incl_b))) + return np.nan_to_num(vt * np.sin(abs(np.radians(self.incl_b)))) # -- The emission cube: returns in [Jy/pixel]; velax and vlsr in [m/s] @@ -250,7 +263,7 @@ def get_cube(self, velax, restfreq=230.538e9, vlsr=0.0): # emission distribution (front surface) Inu_f = Bnu_f * (1.0 - np.exp(-tau_nuf)) - Inu_f = np.where(np.isfinite(Inu_f), Inu_f, 0.0) + Inu_f = np.where(np.isfinite(Inu_f), Inu_f, 0.0) # spectrally-dependent optical depths (back surface) tau_nub = self.gaussian(velax[:, None, None], @@ -265,11 +278,12 @@ def get_cube(self, velax, restfreq=230.538e9, vlsr=0.0): # emission distribution (back surface) Inu_b = Bnu_b * (1.0 - np.exp(-tau_nub)) * np.exp(-tau_nuf) - Inu_b = np.where(np.isfinite(Inu_b), Inu_b, 0.0) + Inu_b = np.where(np.isfinite(Inu_b), Inu_b, 0.0) # combined emission distribution, to proper Jy/pixel units - Inu = Inu_f + Inu_b + Inu = Inu_f + Inu_b pix_area = (self.cell_sky * np.pi / 180. / 3600.)**2 Inu *= 1e26 * pix_area - + + return Inu diff --git a/parametric_disk_CSALT0.py b/parametric_disk_CSALT0.py index 82fbcda..5f0b972 100644 --- a/parametric_disk_CSALT0.py +++ b/parametric_disk_CSALT0.py @@ -8,7 +8,7 @@ def tapered_powerlaw(r, y_0, r_0, y_q, y_tap, y_exp): """Exponentially tapered power law.""" - return y_0 * (r / r_0)**y_q * np.exp(-(r / y_tap)**y_exp) + return np.nan_to_num(y_0 * (r / r_0)**y_q * np.exp(-(r / y_tap)**y_exp)) def parametric_disk(velax, pars, pars_fixed, quiet=True): @@ -36,7 +36,7 @@ def parametric_disk(velax, pars, pars_fixed, quiet=True): dVmax_b = np.sqrt(2 * sc.k * Tmax_b / (28 * (sc.m_p + sc.m_e))) x0, y0 = 0., 0. z_tap, Tb_tap, dV_tap, tau_tap = np.inf, r_l, r_l, r_l - z_exp, Tb_exp, dV_exp, tau_exp = np.inf, np.inf, np.inf, np.inf + z_exp, Tb_exp, dV_exp, tau_exp = np.inf, np.inf, np.inf, np.inf #this functionally sets everything to a power law # Get a simple_disk instance. @@ -85,11 +85,14 @@ def tau(r): # deal with the divergence of power laws close to the disk center. Although # this can be specified for each side independently, we assume they're the # same for simplicity (so setting `side='both'`). - disk.set_Tgas_profile(function=Tgas, min=0.0, max=Tmax_f, side='front') - disk.set_Tgas_profile(function=Tgas, min=0.0, max=Tmax_b, side='back') - disk.set_dV_profile(function=dV, min=0.0, max=dVmax_f, side='front') - disk.set_dV_profile(function=dV, min=0.0, max=dVmax_b, side='back') - disk.set_tau_profile(function=tau, min=0.0, max=None, side='both') + + disk.set_Tgas_profile(function=Tgas, vmin=0.0, vmax=Tmax_f, side='front') + disk.set_Tgas_profile(function=Tgas, vmin=0.0, vmax=Tmax_b, side='back') + disk.set_dV_profile(function=dV, vmin=0.0, vmax=dVmax_f, side='front') + disk.set_dV_profile(function=dV, vmin=0.0, vmax=dVmax_b, side='back') + disk.set_tau_profile(function=tau, vmin=0.0, vmax=None, side='front') + disk.set_tau_profile(function=tau, vmin=0.0, vmax=None, side='back') + # Set up the velocity structure. Here we use a simple Keplerian rotation @@ -99,7 +102,7 @@ def vkep(r): r_m, z_m = r * sc.au, z_f(r / dist) * dist * sc.au vv = np.sqrt(sc.G * mstar * 1.98847e30 * r_m**2 / \ np.power(r_m**2 + z_m**2, 1.5)) - return vv + return np.nan_to_num(vv) disk.set_vtheta_profile(function=vkep, side='both')