diff --git a/vip_hci/fm/scattered_light_disk.py b/vip_hci/fm/scattered_light_disk.py index 7877c37c..f11f9fda 100644 --- a/vip_hci/fm/scattered_light_disk.py +++ b/vip_hci/fm/scattered_light_disk.py @@ -19,7 +19,7 @@ import numpy as np import matplotlib.pyplot as plt from scipy.optimize import newton -from scipy.interpolate import interp1d +from scipy.interpolate import interp1d,PchipInterpolator from ..var import frame_center @@ -36,7 +36,7 @@ def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., 'a': 40, 'e': 0, 'ksi0': 1., 'gamma': 2., 'beta': 1., 'dens_at_r0': 1.}, spf_dico={'name': 'HG', 'g': 0., 'polar': False}, xdo=0., - ydo=0.): + ydo=0.,xs=None,ys=None): """ Constructor of the Scattered_light_disk object, taking in input the geometric parameters of the disk, the radial density distribution @@ -44,7 +44,7 @@ def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., So far, only one radial distribution is implemented: a smoothed 2-power-law distribution, but more complex radial profiles can be implemented on demand. - The star is assumed to be centered at the frame center as defined in + By default, the star is assumed to be centered at the frame center as defined in the vip_hci.var.frame_center function (geometric center of the image, e.g. either in the middle of the central pixel for odd-size images or in between 4 pixel for even-size images). @@ -107,6 +107,18 @@ def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., ydo : float disk offset along the y-axis in the disk frame (=semi-minor axis), in a.u. (default 0) + xs : float + star x position in pixel. By default, the star is assumed to be + centered at the frame center as defined in + the vip_hci.var.frame_center function (geometric center of the image, + e.g. either in the middle of the central pixel for odd-size images or + in between 4 pixel for even-size images). + ys : float + star y position in pixel. By default, the star is assumed to be + centered at the frame center as defined in + the vip_hci.var.frame_center function (geometric center of the image, + e.g. either in the middle of the central pixel for odd-size images or + in between 4 pixel for even-size images). """ self.nx = nx # number of pixels along the x axis of the image self.ny = ny # number of pixels along the y axis of the image @@ -123,7 +135,11 @@ def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., self.rmin = np.sqrt(self.xdo**2+self.ydo**2)+self.pxInAU self.dust_density = Dust_distribution(density_dico) # star center along the y- and x-axis, in pixels - self.yc, self.xc = frame_center(np.ndarray((self.ny, self.nx))) + if xs is None or ys is None: + self.yc, self.xc = frame_center(np.ndarray((self.ny, self.nx))) + else: + self.yc = ys + self.xc = xs self.x_vector = (np.arange(0, nx) - self.xc)*self.pxInAU # x axis in au self.y_vector = (np.arange(0, ny) - self.yc)*self.pxInAU # y axis in au self.x_map_0PA, self.y_map_0PA = np.meshgrid(self.x_vector, @@ -356,6 +372,40 @@ def compute_scattered_light(self, halfNbSlices=25): np.nanmax(self.scattered_light_map)) return self.scattered_light_map + def get_scattering_angle(self): + """ + Function that computes an image with the scattering angle of the dust + located in the disk midplane at each point in the image + + Returns + ------- + numpy.ndarray + image if size nx , ny. + + """ + # dist along the line of sight to reach the disk midplane (z_D=0), AU: + lz0_map = self.y_map * np.tan(np.deg2rad(self.itilt)) + + # 1d array pre-calculated values, AU + ycs_vector = self.cosi*self.y_map + # 1d array pre-calculated values, AU + zsn_vector = -self.sini*self.y_map + xd_vector = self.x_map # x_disk, in AU + # rotation about x axis + yd_vector = ycs_vector + self.sini * lz0_map # y_Disk in AU + zd_vector = zsn_vector + self.cosi * lz0_map # z_Disk, in AU + # Dist and polar angles in the frame centered on the star position: + # squared distance to the star, in AU^2 + d2star_vector = xd_vector**2+yd_vector**2+zd_vector**2 + dstar_vector = np.sqrt(d2star_vector) # distance to the star, in AU + # midplane distance to the star (r coordinate), in AU + rstar_vector = np.sqrt(xd_vector**2+yd_vector**2) + thetastar_vector = np.arctan2(yd_vector, xd_vector) + # Phase angles: + cosphi_vector = (rstar_vector*self.sini*np.sin(thetastar_vector) + + zd_vector*self.cosi)/dstar_vector # in radians + return np.rad2deg(np.arccos(cosphi_vector)) + class Dust_distribution(object): """This class represents the dust distribution @@ -1063,6 +1113,8 @@ def interpolate_phase_function(self, spf_dico): Optionnaly it can specify the kind of interpolation requested (as specified by the scipy.interpolate.interp1d function), by default it uses a quadratic interpolation. + best_spf_interp_pchip = PchipInterpolator(scat_angles,best_spf)(phi) + """ if 'kind' in spf_dico.keys(): if not isinstance(spf_dico['kind'], int) and spf_dico['kind'] not in \ @@ -1071,15 +1123,19 @@ def interpolate_phase_function(self, spf_dico): raise TypeError('The key "{0:s}" must be an integer ' 'or a string ("linear", "nearest", "zero", "slinear", ' '"quadratic", "cubic", "previous",' - '"next"'.format(spf_dico['kind'])) + '"next" or "pchip")'.format(spf_dico['kind'])) else: kind = spf_dico['kind'] else: - kind = 'quadratic' - self.interpolation_function = interp1d(spf_dico['phi'], - spf_dico['spf'], kind=kind, - bounds_error=False, - fill_value=np.nan) + kind = 'pchip' + # by default, we use pchip for interpolation + interp_func = PchipInterpolator(spf_dico['phi'],spf_dico['spf']) + if kind != 'pchip': + interp_func = interp1d(spf_dico['phi'],spf_dico['spf'], kind=kind, + bounds_error=False, + fill_value=np.nan) + + self.interpolation_function = interp_func def compute_phase_function_from_cosphi(self, cos_phi): """