Skip to content

Correcting deprojection problems #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
j6626 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Correcting deprojection problems #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 26 additions & 12 deletions csalt/csalt_disk.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
from scipy.spatial import Delaunay
import scipy.constants as sc
import warnings

from scipy.interpolate import LinearNDInterpolator

class csalt_disk:

Expand Down Expand Up @@ -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)))
Expand All @@ -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,
Expand Down Expand Up @@ -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]
Expand All @@ -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],
Expand All @@ -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
19 changes: 11 additions & 8 deletions parametric_disk_CSALT0.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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')

Expand Down