From 6cee77c154247830f2e6cfe7a6c326be5f20a2de Mon Sep 17 00:00:00 2001 From: jmilou Date: Mon, 21 Feb 2022 22:10:29 +0100 Subject: [PATCH] improvement and bug correction of phase_function.py Improvement: addition of the possibility to model the phase function as a polynomial Bug correction if the SPF is provided as a piecewise function --- vip_hci/metrics/phase_function.py | 51 +++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/vip_hci/metrics/phase_function.py b/vip_hci/metrics/phase_function.py index ca2dd424..2f0c1f98 100644 --- a/vip_hci/metrics/phase_function.py +++ b/vip_hci/metrics/phase_function.py @@ -33,8 +33,12 @@ def __init__(self, spf_dico={'name': 'HG', 'g': 0., 'polar': False}): (single Henyey Greenstein), 'DoubleHG' (double Heyney Greenstein) or 'interpolated' (custom function). The parameter "polar" enables to switch on the polarisation (if set - to True, the default is False). In this case it assumes a Rayleigh - polarised fraction (1-(cos phi)^2) / (1+(cos phi)^2). + to True, the default is False). In this case it assumes either + - a Rayleigh polarised fraction (1-(cos phi)^2) / (1+(cos phi)^2). + if nothing else is specified + - a polynomial if the keyword 'polar_polynom_coeff' is defined + and corresponds to an array of polynomial coefficient, e.g. + [p3,p2,p1,p0] evaluated as np.polyval([p3,p2,p1,p0],np.arange(0, 180, 1)) """ if not isinstance(spf_dico,dict): msg = 'The parameters describing the phase function must be a ' \ @@ -47,12 +51,22 @@ def __init__(self, spf_dico={'name': 'HG', 'g': 0., 'polar': False}): self.type = spf_dico['name'] if 'polar' not in spf_dico.keys(): self.polar = False - elif not isinstance(spf_dico['polar'], bool): - msg = 'The dictionnary describing the polarisation must be a ' \ - 'boolean' - raise TypeError(msg) - else: + else: + if not isinstance(spf_dico['polar'], bool): + msg = 'The dictionnary describing the polarisation must be a ' \ + 'boolean' + raise TypeError(msg) self.polar = spf_dico['polar'] + if 'polar_polynom_coeff' in spf_dico.keys(): + self.polar_polynom = True + if isinstance(spf_dico['polar_polynom_coeff'], (tuple,list,np.ndarray)): + self.polar_polynom_coeff = spf_dico['polar_polynom_coeff'] + else: + msg = 'The dictionnary describing the polarisation polynomial function must be an ' \ + 'array' + raise TypeError(msg) + else: + self.polar_polynom = False if self.type == 'HG': self.phase_function_calc = HenyeyGreenstein_SPF(spf_dico) elif self.type == 'DoubleHG': @@ -77,7 +91,11 @@ def compute_phase_function_from_cosphi(self, cos_phi): """ phf = self.phase_function_calc.compute_phase_function_from_cosphi(cos_phi) if self.polar: - return (1-cos_phi**2)/(1+cos_phi**2) * phf + if self.polar_polynom: + phi = np.rad2deg(np.arccos(cos_phi)) + return np.polyval(self.polar_polynom_coeff, phi) * phf + else: + return (1-cos_phi**2)/(1+cos_phi**2) * phf else: return phf @@ -100,8 +118,11 @@ def plot_phase_function(self): phi = np.arange(0, 180, 1) phase_func = self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi))) if self.polar: - phase_func = (1-np.cos(np.deg2rad(phi))**2) / \ - (1+np.cos(np.deg2rad(phi))**2) * phase_func + if self.polar_polynom: + phase_func = np.polyval(self.polar_polynom_coeff, phi) * phase_func + else: + phase_func = (1-np.cos(np.deg2rad(phi))**2) / \ + (1+np.cos(np.deg2rad(phi))**2) * phase_func plt.close(0) plt.figure(0) @@ -331,10 +352,15 @@ def interpolate_phase_function(self, spf_dico): kind=spf_dico['kind'] else: kind='quadratic' - self.interpolation_function = interp1d(np.cos(np.deg2rad(spf_dico['phi'])),\ + self.interpolation_function = interp1d(spf_dico['phi'],\ spf_dico['spf'],kind=kind,\ bounds_error=False, fill_value=np.nan) + # commented on 2022-02-16 + # self.interpolation_function = interp1d(np.cos(np.deg2rad(spf_dico['phi'])),\ + # spf_dico['spf'],kind=kind,\ + # bounds_error=False, + # fill_value=np.nan) def compute_phase_function_from_cosphi(self, cos_phi): """ @@ -348,4 +374,5 @@ def compute_phase_function_from_cosphi(self, cos_phi): cosine of the scattering angle(s) at which the scattering function must be calculated. """ - return self.interpolation_function(cos_phi) + return self.interpolation_function(np.rad2deg(np.arccos(cos_phi))) + # return self.interpolation_function(cos_phi)