From f397617019429542d819a159fdf91bf5fdf93bbe Mon Sep 17 00:00:00 2001
From: Carlos Gomez <carlosgg33@gmail.com>
Date: Wed, 19 Sep 2018 13:14:12 +0200
Subject: [PATCH] Renaming function create_fake_disk to cube_inject_fakedisk.
 Code cleanup to make it more compliant with PEP8

The following modules were modified: phase_function, scattered_light_disk and dust_distribution modules
---
 vip_hci/metrics/dust_distribution.py    | 201 ++++++++++++---------
 vip_hci/metrics/fakedisk.py             |   5 +-
 vip_hci/metrics/phase_function.py       | 228 ++++++++++++++----------
 vip_hci/metrics/scattered_light_disk.py |   4 +-
 4 files changed, 253 insertions(+), 185 deletions(-)

diff --git a/vip_hci/metrics/dust_distribution.py b/vip_hci/metrics/dust_distribution.py
index 6caf3112..6e56e85a 100644
--- a/vip_hci/metrics/dust_distribution.py
+++ b/vip_hci/metrics/dust_distribution.py
@@ -1,17 +1,18 @@
 # -*- coding: utf-8 -*-
 """
 Created on Wed May  6 17:07:00 2015
-
-@author: jmilli
 """
 import numpy as np
 
-class Dust_distribution:
+_author__ = 'Julien Milli'
+
+
+class Dust_distribution(object):
     """This class represents the dust distribution
     """
-    
-    def __init__(self,density_dico={'name':'2PowerLaws','ain':5,'aout':-5,\
-                               'a':60,'e':0,'ksi0':1.,'gamma':2.,'beta':1.}):
+    def __init__(self,density_dico={'name':'2PowerLaws', 'ain':5, 'aout':-5,
+                                    'a':60, 'e':0, 'ksi0':1., 'gamma':2.,
+                                    'beta':1.}):
         """ 
         Constructor for the Dust_distribution class.
         
@@ -20,55 +21,63 @@ def __init__(self,density_dico={'name':'2PowerLaws','ain':5,'aout':-5,\
         the midplane, and vertically whenever it drops below 0.5% of the 
         peak density in the midplane
         """
-        self.accuracy=5.e-3 
-        if not isinstance(density_dico,dict):
-            raise TypeError('The parameters describing the dust density distribution must be a Python dictionnary')
+        self.accuracy = 5.e-3
+        if not isinstance(density_dico, dict):
+            errmsg = 'The parameters describing the dust density distribution' \
+                     ' must be a Python dictionnary'
+            raise TypeError(errmsg)
         if 'name' not in density_dico.keys():
-            raise TypeError('The dictionnary describing the dust density distribution must contain the key "name"')
+            errmsg = 'The dictionnary describing the dust density ' \
+                     'distribution must contain the key "name"'
+            raise TypeError(errmsg)
         self.type = density_dico['name']
         if self.type == '2PowerLaws':
-            self.dust_distribution_calc = DustEllipticalDistribution2PowerLaws(\
-                                        self.accuracy,density_dico)
+            self.dust_distribution_calc = DustEllipticalDistribution2PowerLaws(
+                                                    self.accuracy, density_dico)
         else:
-            raise TypeError('The only dust distribution implemented so far is the "2PowerLaws"')
+            errmsg = 'The only dust distribution implemented so far is the' \
+                     ' "2PowerLaws"'
+            raise TypeError(errmsg)
 
     def set_density_distribution(self,density_dico):
         """
-        Updates the parameters of the density distribution.
+        Update the parameters of the density distribution.
         """
         self.dust_distribution_calc.set_density_distribution(density_dico)
 
-    def density_cylindrical(self,r,costheta,z):
+    def density_cylindrical(self, r, costheta, z):
         """ 
-        Returns the particule volume density at r, theta, z
+        Return the particule volume density at r, theta, z.
         """
-        return self.dust_distribution_calc.density_cylindrical(r,costheta,z)   
+        return self.dust_distribution_calc.density_cylindrical(r, costheta, z)
 
-    def density_cartesian(self,x,y,z):
+    def density_cartesian(self, x, y, z):
         """ 
-        Returns the particule volume density at x,y,z, taking into account the offset 
-        of the disk
+        Return the particule volume density at x,y,z, taking into account the
+        offset of the disk.
         """           
-        return self.dust_distribution_calc.density_cartesian(x,y,z)
+        return self.dust_distribution_calc.density_cartesian(x, y, z)
 
-        
-    def print_info(self,pxInAu=None):
+    def print_info(self, pxInAu=None):
         """
-        Utility function that displays the parameters of the radial distribution of the dust
+        Utility function that displays the parameters of the radial distribution
+        of the dust
+
         Input:
             - pxInAu (optional): the pixel size in au
         """
         print('----------------------------')
         print('Dust distribution parameters')
         print('----------------------------')
-        self.dust_distribution_calc.print_info(pxInAu=None)
-        
+        self.dust_distribution_calc.print_info(pxInAu)
+
+
 class DustEllipticalDistribution2PowerLaws:
     """
     """
-
-    def __init__(self,accuracy=5.e-3,density_dico={'ain':5,'aout':-5,\
-                           'a':60,'e':0,'ksi0':1.,'gamma':2.,'beta':1.}):
+    def __init__(self, accuracy=5.e-3, density_dico={'ain':5,'aout':-5,
+                                                     'a':60,'e':0,'ksi0':1.,
+                                                     'gamma':2.,'beta':1.}):
         """ 
         Constructor for the Dust_distribution class.
         
@@ -84,37 +93,37 @@ def set_density_distribution(self,density_dico):
         """
         """
         if 'ksi0' not in density_dico.keys():
-            ksi0=1.
+            ksi0 = 1.
         else:
-            ksi0=density_dico['ksi0']
+            ksi0 = density_dico['ksi0']
         if 'beta' not in density_dico.keys():
-            beta=1.
+            beta = 1.
         else:
-            beta=density_dico['beta']
+            beta = density_dico['beta']
         if 'gamma' not in density_dico.keys():
-            gamma=1.
+            gamma = 1.
         else:
-            gamma=density_dico['gamma']
+            gamma = density_dico['gamma']
         if 'aout' not in density_dico.keys():
-            aout=-5.
+            aout = -5.
         else:
-            aout=density_dico['aout']
+            aout = density_dico['aout']
         if 'ain' not in density_dico.keys():
-            ain=5.
+            ain = 5.
         else:
-            ain=density_dico['ain']
+            ain = density_dico['ain']
         if 'e' not in density_dico.keys():
-            e=0.
+            e = 0.
         else:
-            e=density_dico['e']
+            e = density_dico['e']
         if 'a' not in density_dico.keys():
-            a=60.
+            a = 60.
         else:
-            a=density_dico['a']
-        self.set_vertical_density(ksi0=ksi0,gamma=gamma,beta=beta)
-        self.set_radial_density(ain=ain,aout=aout,a=a,e=e)
+            a = density_dico['a']
+        self.set_vertical_density(ksi0=ksi0, gamma=gamma, beta=beta)
+        self.set_radial_density(ain=ain, aout=aout, a=a, e=e)
 
-    def set_vertical_density(self,ksi0=1.,gamma=2.,beta=1.):      
+    def set_vertical_density(self, ksi0=1., gamma=2., beta=1.):
         """ 
         Sets the parameters of the vertical density function    
 
@@ -139,12 +148,12 @@ def set_vertical_density(self,ksi0=1.,gamma=2.,beta=1.):
             print('Warning the flaring coefficient beta is negative')
             print('beta was changed from {0:6.2f} to 0 (flat disk)'.format(beta))
             beta = 0.
-        self.ksi0=float(ksi0)
-        self.gamma=float(gamma)
-        self.beta=float(beta)      
-        self.zmax=ksi0*(-np.log(self.accuracy))**(1./gamma)
+        self.ksi0 = float(ksi0)
+        self.gamma = float(gamma)
+        self.beta = float(beta)
+        self.zmax = ksi0*(-np.log(self.accuracy))**(1./gamma)
 
-    def set_radial_density(self,ain=5.,aout=-5.,a=60.,e=0.):   
+    def set_radial_density(self, ain=5., aout=-5., a=60., e=0.):
         """ 
         Sets the parameters of the radial density function
 
@@ -179,15 +188,17 @@ def set_radial_density(self,ain=5.,aout=-5.,a=60.,e=0.):
             e = 0.99  
         if a < 0:
             raise ValueError('Warning the semi-major axis a is negative')
-        self.ain=float(ain)
-        self.aout=float(aout)
-        self.a=float(a)
-        self.e=float(e)
-        self.p=self.a*(1-self.e**2)
+        self.ain = float(ain)
+        self.aout = float(aout)
+        self.a = float(a)
+        self.e = float(e)
+        self.p = self.a*(1-self.e**2)
         try:
-            self.rmax=self.a*self.accuracy**(1/self.aout) #maximum distance of integration, AU
-            if self.ain!=self.aout:
-                self.apeak = self.a * np.power(-self.ain/self.aout,1./(2.*(self.ain-self.aout)))
+            # maximum distance of integration, AU
+            self.rmax = self.a*self.accuracy**(1/self.aout)
+            if self.ain != self.aout:
+                self.apeak = self.a * np.power(-self.ain/self.aout,
+                                               1./(2.*(self.ain-self.aout)))
             else:
                 self.apeak = self.a
         except OverflowError:
@@ -204,29 +215,41 @@ def set_radial_density(self,ain=5.,aout=-5.,a=60.,e=0.):
             raise ZeroDivisionError
         self.itiltthreshold = np.rad2deg(np.arctan(self.rmax/self.zmax))
                 
-    def print_info(self,pxInAu=None):
+    def print_info(self, pxInAu=None):
         """
-        Utility function that displays the parameters of the radial distribution of the dust
+        Utility function that displays the parameters of the radial distribution
+        of the dust
+
         Input:
             - pxInAu (optional): the pixel size in au
         """
         if pxInAu is not None:
-            print('Reference semi-major axis: {0:.1f}au or {1:.1f}px'.format(self.a,self.a/pxInAu))
-            print('Semi-major axis at maximum dust density: {0:.1f}au or {1:.1f}px (same as ref sma if ain=-aout)'.format(self.apeak,self.apeak/pxInAu))
-            print('Ellipse p parameter: {0:.1f}au or {1:.1f}px'.format(self.p,self.p/pxInAu))
+            msg = 'Reference semi-major axis: {0:.1f}au or {1:.1f}px'
+            print(msg.format(self.a,self.a/pxInAu))
+            msg2 = 'Semi-major axis at maximum dust density: {0:.1f}au or ' \
+                   '{1:.1f}px (same as ref sma if ain=-aout)'
+            print(msg2.format(self.apeak,self.apeak/pxInAu))
+            msg3 = 'Ellipse p parameter: {0:.1f}au or {1:.1f}px'
+            print(msg3.format(self.p,self.p/pxInAu))
         else:
             print('Reference semi-major axis: {0:.1f}au'.format(self.a))
-            print('Semi-major axis at maximum dust density: {0:.1f}au (same as ref sma if ain=-aout)'.format(self.apeak))
+            msg = 'Semi-major axis at maximum dust density: {0:.1f}au (same ' \
+                  'as ref sma if ain=-aout)'
+            print(msg.format(self.apeak))
             print('Ellipse p parameter: {0:.1f}au'.format(self.p))
         print('Ellipticity: {0:.3f}'.format(self.e))
         print('Inner slope: {0:.2f}'.format(self.ain))
         print('Outer slope: {0:.2f}'.format(self.aout))
         if pxInAu is not None:
-            print('Scale height: {0:.1f}au or {1:.1f}px at {2:.1f}'.format(self.ksi0,self.ksi0/pxInAu,self.a))
+            msg = 'Scale height: {0:.1f}au or {1:.1f}px at {2:.1f}'
+            print(msg.format(self.ksi0,self.ksi0/pxInAu,self.a))
         else:
-            print('Scale height: {0:.2f} au at {1:.2f}'.format(self.ksi0,self.a))            
+            print('Scale height: {0:.2f} au at {1:.2f}'.format(self.ksi0,
+                                                               self.a))
         print('Vertical profile index: {0:.2f}'.format(self.gamma))
-        print('Disc vertical FWHM: {0:.2f} at {1:.2f}'.format(2.*self.ksi0*np.power(np.log10(2.),1./self.gamma),self.a))
+        msg = 'Disc vertical FWHM: {0:.2f} at {1:.2f}'
+        print(msg.format(2.*self.ksi0*np.power(np.log10(2.), 1./self.gamma),
+                         self.a))
         print('Flaring coefficient: {0:.2f}'.format(self.beta))
         print('------------------------------------')
         print('Properties for numerical integration')
@@ -234,39 +257,43 @@ def print_info(self,pxInAu=None):
         print('Requested accuracy {0:.2e}'.format(self.accuracy))
 #        print('Minimum radius for integration: {0:.2f} au'.format(self.rmin))
         print('Maximum radius for integration: {0:.2f} au'.format(self.rmax))
-        print('Maximum height for integration: {0:.2f} au'.format(self.zmax))      
-        print('Inclination threshold: {0:.2f} degrees'.format(self.itiltthreshold))
+        print('Maximum height for integration: {0:.2f} au'.format(self.zmax))
+        msg = 'Inclination threshold: {0:.2f} degrees'
+        print(msg.format(self.itiltthreshold))
         return
         
-    def density_cylindrical(self,r,costheta,z):
+    def density_cylindrical(self, r, costheta, z):
         """ Returns the particule volume density at r, theta, z
         """
-        radial_ratio=r/(self.p/(1-self.e*costheta))
-        radial_density_term=np.sqrt(2./(np.power(radial_ratio,-2*self.ain)+np.power(radial_ratio,-2*self.aout)))
-        vertical_density_term=np.exp(-np.power(np.abs(z)/(self.ksi0*np.power(radial_ratio,self.beta)),self.gamma))
+        radial_ratio = r/(self.p/(1-self.e*costheta))
+        den = (np.power(radial_ratio, -2*self.ain) +
+               np.power(radial_ratio,-2*self.aout))
+        radial_density_term = np.sqrt(2./den)
+        den2 = (self.ksi0*np.power(radial_ratio,self.beta))
+        vertical_density_term = np.exp(-np.power(np.abs(z)/den2, self.gamma))
         return radial_density_term*vertical_density_term    
 
-    def density_cartesian(self,x,y,z):
-        """ Returns the particule volume density at x,y,z, taking into account the offset 
-        of the disk
+    def density_cartesian(self, x, y, z):
+        """ Returns the particule volume density at x,y,z, taking into account
+        the offset of the disk
         """
-        r=np.sqrt(x**2+y**2)
-        if r==0:
+        r = np.sqrt(x**2+y**2)
+        if r == 0:
             costheta=0
         else:
-            costheta=x/r            
+            costheta = x/r
         return self.density_cylindrical(r,costheta,z)
 
-if __name__=='__main__':
+
+if __name__ == '__main__':
     """
     Main of the class for debugging
     """
     test = DustEllipticalDistribution2PowerLaws()
-    test.set_radial_density(ain=5.,aout=-5.,a=60.,e=0.)
+    test.set_radial_density(ain=5., aout=-5., a=60., e=0.)
     test.print_info()
-    costheta=0.
-    z=0.
+    costheta = 0.
+    z = 0.
     for a in np.linspace(60-5,60+5,11):
-        t = test.density_cylindrical(a,costheta,z)
-        print('r={0:.1f} density={1:.4f}'.format(a,t))
-    
\ No newline at end of file
+        t = test.density_cylindrical(a, costheta, z)
+        print('r={0:.1f} density={1:.4f}'.format(a, t))
\ No newline at end of file
diff --git a/vip_hci/metrics/fakedisk.py b/vip_hci/metrics/fakedisk.py
index dd5d44ce..9c6714dd 100644
--- a/vip_hci/metrics/fakedisk.py
+++ b/vip_hci/metrics/fakedisk.py
@@ -7,17 +7,16 @@
 from __future__ import division, print_function
 
 __author__ = 'Julien Milli @ ESO, Valentin Christiaens @ ULg/UChile'
-__all__ = ['create_fakedisk_cube',
+__all__ = ['cube_inject_fakedisk',
            'cube_inject_trace']
 
 import numpy as np
 from scipy import signal
 from ..preproc import cube_derotate, frame_shift
 from ..var import frame_center
-import pdb
 
 
-def create_fakedisk_cube(fakedisk, angle_list, psf=None, imlib='opencv',
+def cube_inject_fakedisk(fakedisk, angle_list, psf=None, imlib='opencv',
                          interpolation='lanczos4', cxy=None, nproc=1,
                          border_mode='constant'):
     """
diff --git a/vip_hci/metrics/phase_function.py b/vip_hci/metrics/phase_function.py
index da681148..d9f99a51 100644
--- a/vip_hci/metrics/phase_function.py
+++ b/vip_hci/metrics/phase_function.py
@@ -6,21 +6,22 @@
 from __future__ import division, print_function
 
 __author__ = 'Julien Milli'
-__all__ = [] #['create_fakedisk_cube']
+__all__ = []
 
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.interpolate import interp1d
 
-class Phase_function:
-    """This class represents the scattering phase function (spf). 
+
+class Phase_function(object):
+    """ This class represents the scattering phase function (spf).
     It contains the attribute phase_function_calc that implements either a 
     single Henyey Greenstein phase function, a double Heyney Greenstein, 
     or any custom function (data interpolated from 
     an input list of phi, spf(phi)).
     """
     
-    def __init__(self,spf_dico={'name':'HG','g':0.,'polar':False}):
+    def __init__(self, spf_dico={'name': 'HG', 'g': 0., 'polar': False}):
         """
         Constructor of the Phase_function class. It checks whether the spf_dico
         contains a correct name and sets the attribute phase_function_calc
@@ -28,24 +29,30 @@ def __init__(self,spf_dico={'name':'HG','g':0.,'polar':False}):
         Parameters
         ----------
         spf_dico :  dictionnary
-            Parameters describing the scattering phase function to be implemented.
-            By default, an isotropic phase function is implemented. 
-            It should at least contain the key "name" chosen between 'HG' (single Henyey 
-            Greenstein), 'DoubleHG' (double Heyney Greenstein) or
+            Parameters describing the scattering phase function to be
+            implemented. By default, an isotropic phase function is implemented.
+            It should at least contain the key "name" chosen between 'HG'
+            (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)            
+            polarised fraction (1-(cos phi)^2) / (1+(cos phi)^2).
         """
         if not isinstance(spf_dico,dict):
-            raise TypeError('The parameters describing the phase function must be a Python dictionnary')
+            msg = 'The parameters describing the phase function must be a ' \
+                  'Python dictionnary'
+            raise TypeError(msg)
         if 'name' not in spf_dico.keys():
-            raise TypeError('The dictionnary describing the phase function must contain the key "name"')
+            msg = 'The dictionnary describing the phase function must contain' \
+                  ' the key "name"'
+            raise TypeError(msg)
         self.type = spf_dico['name']
         if 'polar' not in spf_dico.keys():
-            self.polar=False
-        elif not isinstance(spf_dico['polar'],bool):
-            raise TypeError('The dictionnary describing the polarisation must be a boolean')
+            self.polar = False
+        elif not isinstance(spf_dico['polar'], bool):
+            msg = 'The dictionnary describing the polarisation must be a ' \
+                  'boolean'
+            raise TypeError(msg)
         else: 
             self.polar = spf_dico['polar']
         if self.type == 'HG':
@@ -55,24 +62,26 @@ def __init__(self,spf_dico={'name':'HG','g':0.,'polar':False}):
         elif self.type == 'interpolated':
             self.phase_function_calc = Interpolated_SPF(spf_dico)
         else:
-            raise TypeError('Type of phase function not understood: {0:s}'.format(self.type))
+            msg = 'Type of phase function not understood: {0:s}'
+            raise TypeError(msg.format(self.type))
                 
-    def compute_phase_function_from_cosphi(self,cos_phi):
+    def compute_phase_function_from_cosphi(self, cos_phi):
         """
-        Compute the phase function at (a) specific scattering scattering angle(s) 
-        phi. The argument is not phi but cos(phi) for optimization reasons.
+        Compute the phase function at (a) specific scattering scattering
+        angle(s) phi. The argument is not phi but cos(phi) for optimization
+        reasons.
 
         Parameters
         ----------
         cos_phi : float or array
-            cosine of the scattering angle(s) at which the scattering function must
-            be calculated.
+            cosine of the scattering angle(s) at which the scattering function
+            must be calculated.
         """
+        phf = self.phase_function_calc.compute_phase_function_from_cosphi(cos_phi)
         if self.polar:
-            return (1-cos_phi**2)/(1+cos_phi**2) * \
-                self.phase_function_calc.compute_phase_function_from_cosphi(cos_phi)
+            return (1-cos_phi**2)/(1+cos_phi**2) * phf
         else:
-            return self.phase_function_calc.compute_phase_function_from_cosphi(cos_phi)
+            return phf
 
     def print_info(self):
         """
@@ -90,28 +99,29 @@ def plot_phase_function(self):
         """
         Plots the scattering phase function
         """
-        phi = np.arange(0,180,1)
+        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) * \
-                self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi)))
-        else:
-            phase_func = self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi)))
+            phase_func = (1-np.cos(np.deg2rad(phi))**2) / \
+                         (1+np.cos(np.deg2rad(phi))**2) * phase_func
+
         plt.close(0)
         plt.figure(0)
-        plt.plot(phi,phase_func)
+        plt.plot(phi, phase_func)
         plt.xlabel('Scattering phase angle in degrees')
         plt.ylabel('Scattering phase function')
         plt.grid()
-        plt.xlim(0,180)
+        plt.xlim(0, 180)
         plt.show()
 
-class HenyeyGreenstein_SPF:
+
+class HenyeyGreenstein_SPF(object):
     """
-    Implementation of a scattering phase function with a single Henyey Greenstein 
-    function. 
+    Implementation of a scattering phase function with a single Henyey
+    Greenstein function.
     """
 
-    def __init__(self,spf_dico={'g':0.}):
+    def __init__(self, spf_dico={'g':0.}):
         """
         Constructor of a Heyney Greenstein phase function.
     
@@ -123,36 +133,41 @@ def __init__(self,spf_dico={'g':0.}):
         """
         # it must contain the key "g"
         if 'g' not in spf_dico.keys():
-            raise TypeError('The dictionnary describing a Heyney Greenstein phase function must contain the key "g"')
+            raise TypeError('The dictionnary describing a Heyney Greenstein '
+                            'phase function must contain the key "g"')
         # the value of "g" must be a float or a list of floats
-        elif not isinstance(spf_dico['g'],(float,int)):
-            raise TypeError('The key "g" of a Heyney Greenstein phase function dictionnary must be a float or an integer')
+        elif not isinstance(spf_dico['g'], (float, int)):
+            raise TypeError('The key "g" of a Heyney Greenstein phase function'
+                            ' dictionnary must be a float or an integer')
         self.set_phase_function(spf_dico['g'])
 
-    def set_phase_function(self,g):
+    def set_phase_function(self, g):
         """ 
         Set the value of g
         """
         if g >= 1:
-            print('Warning the Henyey Greenstein parameter is greater than or equal to 1')
+            print('Warning the Henyey Greenstein parameter is greater than or '
+                  'equal to 1')
             print('The value was changed from {0:6.2f} to 0.99'.format(g))
-            g=0.99
+            g = 0.99
         elif g <= -1:
-            print('Warning the Henyey Greenstein parameter is smaller than or equal to -1')
+            print('Warning the Henyey Greenstein parameter is smaller than or '
+                  'equal to -1')
             print('The value was changed from {0:6.2f} to -0.99'.format(g))
-            g=-0.99            
-        self.g=float(g)
+            g = -0.99
+        self.g = float(g)
         
-    def compute_phase_function_from_cosphi(self,cos_phi):
+    def compute_phase_function_from_cosphi(self, cos_phi):
         """
-        Compute the phase function at (a) specific scattering scattering angle(s) 
-        phi. The argument is not phi but cos(phi) for optimization reasons.
+        Compute the phase function at (a) specific scattering scattering
+        angle(s) phi. The argument is not phi but cos(phi) for optimization
+        reasons.
 
         Parameters
         ----------
         cos_phi : float or array
-            cosine of the scattering angle(s) at which the scattering function must
-            be calculated.
+            cosine of the scattering angle(s) at which the scattering function
+            must be calculated.
         """
         return 1./(4*np.pi)*(1-self.g**2)/(1+self.g**2-2*self.g*cos_phi)**(3./2.)         
 
@@ -162,31 +177,42 @@ def print_info(self):
         """
         print('Heynyey Greenstein coefficient: {0:.2f}'.format(self.g))
 
-class DoubleHenyeyGreenstein_SPF:
+
+class DoubleHenyeyGreenstein_SPF(object):
     """
-    Implementation of a scattering phase function with a double Henyey Greenstein 
-    function. 
+    Implementation of a scattering phase function with a double Henyey
+    Greenstein function.
     """
     
-    def __init__(self,spf_dico={'g':[0.5,-0.3],'weight':0.7}):
+    def __init__(self, spf_dico={'g': [0.5,-0.3], 'weight': 0.7}):
         """
         """        
         # it must contain the key "g"
         if 'g' not in spf_dico.keys():
-            raise TypeError('The dictionnary describing a Heyney Greenstein phase function must contain the key "g"')
+            raise TypeError('The dictionnary describing a Heyney Greenstein'
+                            ' phase function must contain the key "g"')
         # the value of "g" must be a list of floats
         elif not isinstance(spf_dico['g'],(list,tuple,np.ndarray)):
-            raise TypeError('The key "g" of a Heyney Greenstein phase function dictionnary must be  a list of floats')
+            raise TypeError('The key "g" of a Heyney Greenstein phase '
+                            'function dictionnary must be  a list of floats')
         # it must contain the key "weight"
         if 'weight' not in spf_dico.keys():
-                raise TypeError('The dictionnary describing a multiple Henyey Greenstein phase function must contain the key "weight"')
+                raise TypeError('The dictionnary describing a multiple Henyey '
+                                'Greenstein phase function must contain the '
+                                'key "weight"')
         # the value of "weight" must be a list of floats
-        elif not isinstance(spf_dico['weight'],(float,int)):
-            raise TypeError('The key "weight" of a Heyney Greenstein phase function dictionnary must be a float (weight of the first HG coefficient between 0 and 1)')
+        elif not isinstance(spf_dico['weight'], (float, int)):
+            raise TypeError('The key "weight" of a Heyney Greenstein phase '
+                            'function dictionnary must be a float (weight of '
+                            'the first HG coefficient between 0 and 1)')
         elif spf_dico['weight']<0 or spf_dico['weight']>1:
-            raise ValueError('The key "weight" of a Heyney Greenstein phase function dictionnary must be between 0 and 1. It corresponds to the weight of the first HG coefficient')            
+            raise ValueError('The key "weight" of a Heyney Greenstein phase'
+                             ' function dictionnary must be between 0 and 1. It'
+                             ' corresponds to the weight of the first HG '
+                             'coefficient')
         if len(spf_dico['g']) != 2:
-            raise TypeError('The keys "weight" and "g" must contain the same number of elements')
+            raise TypeError('The keys "weight" and "g" must contain the same'
+                            ' number of elements')
         self.g = spf_dico['g']
         self.weight = spf_dico['weight']   
 
@@ -194,48 +220,56 @@ def print_info(self):
         """
         Prints the value of the HG coefficients and weights
         """
-        print('Heynyey Greenstein first component : coeff {0:.2f} , weight {1:.1f}%'.format(self.g[0],self.weight*100))
-        print('Heynyey Greenstein second component: coeff {0:.2f} , weight {1:.1f}%'.format(self.g[1],(1-self.weight)*100.))
+        print('Heynyey Greenstein first component : coeff {0:.2f} , '
+              'weight {1:.1f}%'.format(self.g[0], self.weight*100))
+        print('Heynyey Greenstein second component: coeff {0:.2f} , '
+              'weight {1:.1f}%'.format(self.g[1], (1-self.weight)*100.))
 
-    def compute_singleHG_from_cosphi(self,g,cos_phi):
+    def compute_singleHG_from_cosphi(self, g, cos_phi):
         """
-        Compute a single Heyney Greenstein phase function at (a) specific scattering scattering angle(s) 
-        phi. The argument is not phi but cos(phi) for optimization reasons.
+        Compute a single Heyney Greenstein phase function at (a) specific
+        scattering scattering angle(s) phi. The argument is not phi but cos(phi)
+        for optimization reasons.
 
         Parameters
         ----------
         g : float
             Heyney Greenstein coefficient
         cos_phi : float or array
-            cosine of the scattering angle(s) at which the scattering function must
-            be calculated.
+            cosine of the scattering angle(s) at which the scattering function
+            must be calculated.
         """
         return 1./(4*np.pi)*(1-g**2)/(1+g**2-2*g*cos_phi)**(3./2.)         
 
     def compute_phase_function_from_cosphi(self,cos_phi):
         """
-        Compute the phase function at (a) specific scattering scattering angle(s) 
-        phi. The argument is not phi but cos(phi) for optimization reasons.
+        Compute the phase function at (a) specific scattering scattering
+        angle(s) phi. The argument is not phi but cos(phi) for optimization
+        reasons.
 
         Parameters
         ----------
         cos_phi : float or array
-            cosine of the scattering angle(s) at which the scattering function must
-            be calculated.
+            cosine of the scattering angle(s) at which the scattering function
+            must be calculated.
         """
-        return self.weight * self.compute_singleHG_from_cosphi(self.g[0],cos_phi) + \
-                (1-self.weight) * self.compute_singleHG_from_cosphi(self.g[1],cos_phi)
+        return self.weight * self.compute_singleHG_from_cosphi(self.g[0],
+                                                               cos_phi) + \
+                (1-self.weight) * self.compute_singleHG_from_cosphi(self.g[1],
+                                                                    cos_phi)
 
-class Interpolated_SPF:
+
+class Interpolated_SPF(object):
     """
     Custom implementation of a scattering phase function by providing a list of 
     scattering phase angles and corresponding values of the phase function. 
     """
     
-    def __init__(self,spf_dico={'phi':\
-                np.array([  0,  18,  36,  54,  72,  90, 108, 126, 144, 162]),\
-                'spf':np.array([3.580, 0.703, 0.141, 0.0489, 0.0233,\
-                                0.0136, 0.0091, 0.0069, 0.0056,0.005])}):
+    def __init__(self, spf_dico={'phi':np.array([  0,  18,  36,  54,  72,  90,
+                                                   108, 126, 144, 162]),
+                                 'spf':np.array([3.580, 0.703, 0.141, 0.0489,
+                                                 0.0233, 0.0136, 0.0091, 0.0069,
+                                                 0.0056,0.005])}):
         """
         Constructor of the Interpolated_SPF class. It checks whether the spf_dico
         contains the keys 'phi' and 'spf'
@@ -248,26 +282,32 @@ def __init__(self,spf_dico={'phi':\
         """
         for key in ['phi','spf']:
             if key not in spf_dico.keys():
-                raise TypeError('The dictionnary describing a "interpolated" phase function must contain the key "{0:s}"'.format(key))
+                raise TypeError('The dictionnary describing a '
+                                '"interpolated" phase function must contain '
+                                'the key "{0:s}"'.format(key))
             elif isinstance(spf_dico[key],(list,tuple,np.ndarray)):
-                raise TypeError('The key "{0:s}" of a "interpolated" phase function dictionnary must be a list, array or tuple'.format(key))
+                raise TypeError('The key "{0:s}" of a "interpolated" phase'
+                                ' function dictionnary must be a list, array'
+                                ' or tuple'.format(key))
         if len(spf_dico['phi']) != len(spf_dico['spf']):
-            raise TypeError('The keys "phi" and "spf" must contain the same number of elements')
+            raise TypeError('The keys "phi" and "spf" must contain the same'
+                            ' number of elements')
         self.interpolate_phase_function(spf_dico)
 
     def print_info(self):
         """
         Prints the information of the spf
         """
-        phi = np.linspace(0,180,19)
+        phi = np.linspace(0, 180, 19)
         spf = self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi)))
-        print('Scattering angle: ',phi)
-        print('Interpolated scattering phase function: ',spf)
+        print('Scattering angle: ', phi)
+        print('Interpolated scattering phase function: ', spf)
 
-    def interpolate_phase_function(self,spf_dico):
+    def interpolate_phase_function(self, spf_dico):
         """
         Creates the function that returns the scattering phase function based
-        on the scattering angle by interpolating the values given in the dictionnary.
+        on the scattering angle by interpolating the values given in the
+        dictionnary.
     
         Parameters
         ----------
@@ -275,18 +315,20 @@ def interpolate_phase_function(self,spf_dico):
             dictionnary containing the keys "phi" (list of scattering angles)
             and "spf" (list of corresponding scattering phase function values)
         """
-        self.interpolation_function = interp1d(spf_dico['phi'],spf_dico['spf'],\
-                        kind='cubic',bounds_error=False,fill_value=np.nan)
+        self.interpolation_function = interp1d(spf_dico['phi'], spf_dico['spf'],
+                                               kind='cubic', bounds_error=False,
+                                               fill_value=np.nan)
 
-    def compute_phase_function_from_cosphi(self,cos_phi):
+    def compute_phase_function_from_cosphi(self, cos_phi):
         """
-        Compute the phase function at (a) specific scattering scattering angle(s) 
-        phi. The argument is not phi but cos(phi) for optimization reasons.
+        Compute the phase function at (a) specific scattering scattering
+        angle(s) phi. The argument is not phi but cos(phi) for optimization
+        reasons.
 
         Parameters
         ----------
         cos_phi : float or array
-            cosine of the scattering angle(s) at which the scattering function must
-            be calculated.
+            cosine of the scattering angle(s) at which the scattering function
+            must be calculated.
         """
-        return self.interpolation_function(cos_phi)
\ No newline at end of file
+        return self.interpolation_function(cos_phi)
diff --git a/vip_hci/metrics/scattered_light_disk.py b/vip_hci/metrics/scattered_light_disk.py
index 0ff6a91e..280d43df 100644
--- a/vip_hci/metrics/scattered_light_disk.py
+++ b/vip_hci/metrics/scattered_light_disk.py
@@ -323,8 +323,8 @@ def compute_scattered_light(self, halfNbSlices=25):
             limage[il,:,:] = image
         self.scattered_light_map.fill(0.)
         for il in range(1,nbSlices):
-            self.scattered_light_map += (ll[il]-ll[il-1]) * (limage[il-1,:,:] \
-                                                             +limage[il,:,:])
+            self.scattered_light_map += (ll[il]-ll[il-1]) * (limage[il-1,:,:] +
+                                                             limage[il,:,:])
         self.scattered_light_map *= (self.flux_max/np.nanmax(self.scattered_light_map))
         return self.scattered_light_map