<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 6. ADI forward modeling of disks
Author
Step1:  In the following box we import all the VIP routines that will be used in this tutorial.
The path to some routines has changed between versions 1.0.3 and 1.1.0, which saw a major revamp of the modular architecture, hence the if statements.
Step2:  6.1. Introduction
6.1.1. Overview
The functions implemented in vip_hci for disks are located in vip.metrics.scattered_light_disk. It contains the definition of a class called ScatteredLightDisk which can produce a synthetic image of a disk, and also utility functions to create cubes of images where a synthetic disk has been injected at specific position angles to simulate a real observation. 
Currently there is no utility function to do forward modelling and try to find the best disk matching a given dataset as this is usually specific to each dataset.
Keep in mind that ScatteredLightDisk is only a ray-tracing approach and does not contain any physics in it (no radiative transfer, no particle cross-section). It assumes the particles number density around a star follows the mathematical prescription given in section 1.2 and uses a unity scattering cross-section for all particles (no particle size distribution and cross-section dependant on the particle size, the flux of the synthetic disk cannot be converted in physical units (e.g. Jy)
6.1.2. Parametrisation of the density distribution of dust
The density distribution of dust particles is parametrized in a cylindrical coordinate system $\rho(r,\theta,z)$ and is described by the equation
Step3:  6.2.1. Symmetric pole-on disk
For a pole-on disk, $i_\text{tilt}=0^\circ$.
For a symmetric disk, $e=0$ and the position angle (pa) and argument of pericenter ($\omega$) have no impact.
We choose a semi-major axis of 70 a.u., a vertical profile with a gaussian distribution ($\gamma=2$), a reference scale height of 3 a.u. at the semi-major axis of the disk, and inner and outer exponent $\alpha_{in}=12$ and $\alpha_{out}=-12$
Step4:  Then create your disk model
Step5:  The method compute_scattered_light returns the synthetic image of the disk.
Step6:  You can print some info on the geometrical properties of the model, the dust distribution parameters, the numerical integration parameters and the phase function parameters (detailed later). 
This can be useful because, in addition to reminding all the parameters used in the model, it also computes some properties such as the radial FWHM of the disk.
Step7:  As a side note, if $\alpha_{in} \ne \alpha_{out}$, then the peak surface density of the disk is not located at the reference radius $a$.
Step8:  6.2.2. Inclined symmetric disk
Step9:  The position angle of the disk is 0 (e.g. north).  The phase function is asymmetric, the reason why the north and south ansae appear brighter is because the disk is not flat
Step10:  Warning ! The code does not handle perfectly edge-on disks. There is a maximum inclination close to edge-on beyond which it cannot create an image. In practice this is not a limitation as the convolution by the PSF always makes it impossible to disentangle between a close to edge-on disk and a perfectly edge-on disk.
Step11:  6.2.3. Inclined symmetric disk with anisotropy of scattering
6.2.3.1. Simple Henyey-Greenstein phase function
We parametrize the phase function by a Henyey Greenstein phase function, with an asymmetry parameter g. An isotropic phase function has $g=0$, forward scattering is represented by $0<g\leq1$ and backward scattering is represented by $-1\leq g<0$
Step12:  You can plot how the phase function look like
Step13:  The forward side is brighter. 
6.2.3.2. Double Henyey-Greenstein phase function
A double Henyey Greenstein (HG) phase function is simply a linear combination of 2 simple HG phase function. It is therefore parametrized by $g_1$ and $g_2$ the 2 asymmetry parameters of each HG and the weight (between 0 and 1) of the first HG phase function. Typically a double HG is used to represent a combination of forward scattering ($g_1>0$) and backward scattering ($g_2<1$)
Step14:  6.2.3.3. Custom phase function
In some cases, a HG phase function (simple or double) cannot represent well the behaviour of the dust. The code is modular and you can propose new prescriptions for the phase functions if you need, or you can also create a custom phase function.
Step15:  6.2.3.4. Representing a polarised phase function
If you are trying to reproduce the polarised intensity of a disk (for instance Stokes $Q_\phi$ image), you may want to add on top of the scattering phase function, a modulation representing the degree of linear polarisation.
This can be done by setting the polar keyword to True and in this case, the model assumes a Rayleigh-like degree of linear polarisation parametrized by $(1-(\cos \phi)^2) / (1+(\cos \phi)^2)$ where $\phi$ is the scattering angle.
Step16:  You can combine this Rayleigh-like degree of linear polarisation with any phase function (simple HG, double HG or custom type).
6.2.4. Asymmetric disk
Be careful here !
There is no consensus in the community on how to parametrize an eccentric dust distribution, so keep in mind that the convention described in section 1.2 is only one way to do so, but does not mean the dust density distribution in an eccentric disk follows this prescription. For instance, around the pericenter particle velocities are higher and one expects more collision to happen which can create an overdensity of particles compared to other regions of the disk. Conversely, particles stay longer at the apocenter because of Kepler's third law, which means that one could also expect a higher density at apocenter... All these physical phenomena are not described in this model.
Let's start woth a pole-on disk to be insensitive to phase function effects
Step17:  The brightness asymmetry here is entirely due to the fact that the brightness at one point in the disk is inversely proportional to the squared distance to the star.
Once you incline the disk, you start seeing the competing effect of the phase function and eccentricity.
Step18:  6.3. Forward modeling of disks
Let's start from our inclined simple HG symmeric disk fake_disk3_map and assume we observe this disk as part of an ADI sequence of 30 images
Step19:  cube_fake_disk3 is now a cube of 30 frames, where the disk has been injected at the correct position angle.
Step20:  Let's visualize the first, middle and last image of the cube.
Step21:  We can now process this cube with median-ADI for instance
Step22:  The example above shows a typical bias that can be induced by ADI on extended disk signals (Milli et al. 2012).
So far we have not dealt with convolution effects. In practice the image of a disk is convolved by the instrumental PSF. 
Let's assume here an instrument having a gaussian PSF with FWHM = 4px, and create a synthetic PSF using the create_synth_psf function
Step23:  Then we inject the disk in the cube and convolve each frame by the PSF
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from hciplot import plot_frames, plot_cubes
from matplotlib.pyplot import *
from matplotlib import pyplot as plt
import numpy as np
from packaging import version
Explanation: 6. ADI forward modeling of disks
Author: Julien Milli
Last update: 23/03/2022
Suitable for VIP v1.0.0 onwards.
Table of contents
6.1. Introduction
6.1.1. Overview
6.1.2. Parametrisation of the density distribution of dust
6.2. Examples of disks
6.2.1. Symmetric pole-on disk
6.2.2. Inclined symmetric disk
6.2.3. Inclined symmetric disk with anisotropy of scattering
6.2.3.1. Simple Henyey-Greenstein phase function
6.2.3.2. Double Henyey-Greenstein phase function
6.2.3.3. Custom phase function    
6.2.3.4. Representing a polarised phase function    
6.2.4. Asymmetric disk
6.3. Forward modeling of disks
This tutorial shows:
how to generate different models of synthetic (debris) disks;
how to inject model disks in ADI cubes, for forward modeling.
Let's first import a couple of external packages needed in this tutorial:
End of explanation
import vip_hci as vip
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) < version.parse("1.0.0"):
    msg = "Please upgrade your version of VIP"
    msg+= "It should be 1.0.0 or above to run this notebook."
    raise ValueError(msg)
elif version.parse(vvip) <= version.parse("1.0.3"):
    from vip_hci.conf import time_ini, timing
    from vip_hci.medsub import median_sub
    from vip_hci.metrics import cube_inject_fakedisk, ScatteredLightDisk
else:
    from vip_hci.config import time_ini, timing
    from vip_hci.fm import cube_inject_fakedisk, ScatteredLightDisk
    from vip_hci.psfsub import median_sub
# common to all versions:
from vip_hci.var import create_synth_psf
Explanation: In the following box we import all the VIP routines that will be used in this tutorial.
The path to some routines has changed between versions 1.0.3 and 1.1.0, which saw a major revamp of the modular architecture, hence the if statements.
End of explanation
pixel_scale=0.01225 # pixel scale in arcsec/px
dstar= 80 # distance to the star in pc
nx = 200 # number of pixels of your image in X
ny = 200 # number of pixels of your image in Y
Explanation: 6.1. Introduction
6.1.1. Overview
The functions implemented in vip_hci for disks are located in vip.metrics.scattered_light_disk. It contains the definition of a class called ScatteredLightDisk which can produce a synthetic image of a disk, and also utility functions to create cubes of images where a synthetic disk has been injected at specific position angles to simulate a real observation. 
Currently there is no utility function to do forward modelling and try to find the best disk matching a given dataset as this is usually specific to each dataset.
Keep in mind that ScatteredLightDisk is only a ray-tracing approach and does not contain any physics in it (no radiative transfer, no particle cross-section). It assumes the particles number density around a star follows the mathematical prescription given in section 1.2 and uses a unity scattering cross-section for all particles (no particle size distribution and cross-section dependant on the particle size, the flux of the synthetic disk cannot be converted in physical units (e.g. Jy)
6.1.2. Parametrisation of the density distribution of dust
The density distribution of dust particles is parametrized in a cylindrical coordinate system $\rho(r,\theta,z)$ and is described by the equation:
$\rho(r,\theta,z) = \rho_0 \times \left( \frac{2}{\left( \frac{r}{R(\theta)} \right)^{-2a_{in}} + \left( \frac{r}{R(\theta)} \right)^{-2a_{out}} }\right)^{1/2} \times e^{\left[ -\left( \frac{z}{H(r) }\right)^\gamma \right]}$
where $R(\theta)$ is called the reference radius. It is simply the radius of the disk $a$ if the dust distribution is centrally symmetric (no eccentricity). If the disk is eccentric, then $R(\theta)$ depends on $\theta$ and is given by the equation of an ellipse in polar coordinates: $R(\theta) = \frac{a(1-e^2)}{1+e \cos{\theta}}$
This equation for $\rho(r,\theta,z)$ is the product of 3 terms:
1. a constant $\rho_0$ which is the surfacce density of the dust in the midplane, at the reference radius $R(\theta)$. 
2. the density distribution in the midplane $z=0$ defined as $\left( \frac{2}{\left( \frac{r}{R(\theta)} \right)^{-2a_{in}} + \left( \frac{r}{R(\theta)} \right)^{-2a_{out}} }\right)^{1/2}$. Such a function ensures that when $r\ll R(\theta)$ then the term is $\propto r^{\alpha_{in}}$ (and we typically use $\alpha_{in}>0$) and when $r\gg R(\theta)$ then the term is $\propto r^{\alpha_{out}}$ (and we typically use $\alpha_{out}<0$).
3. the vertical profile $e^{\left[ -\left( \frac{z}{H(r) }\right)^\gamma \right]}$ is parametrized by an exponential decay of exponent $\gamma$ and scale height $H(r)$. If $\gamma=2$, the vertical profile is Gaussian (and $H(r)$ is proportional to the $\sigma$ or FWHM of the Gaussian (but not strictly equal to any of them). The scale height is further defined as $H(r)=\xi_0 \times \left( \frac{r}{R(\theta)} \right)^\beta$ where $\xi_0$ is the reference scale height at the reference radius $R(\theta)$ and $\beta$ is the flaring coeffient ($\beta=1$ means a linear flaring: the scale height increases linearly with radius).
6.2. Examples of disks
Let's assume we want to create a synthetic image of 200px, containing a disk around a star located at 80 a.u., observed with SPHERE/IRDIS (pixel scale 12.25 mas).
End of explanation
itilt = 0. # inclination of your disk in degrees
a = 70. # semimajoraxis of the disk in au 
ksi0 = 3. # reference scale height at the semi-major axis of the disk
gamma = 2. # exponant of the vertical exponential decay
alpha_in = 12
alpha_out = -12
beta = 1
Explanation: 6.2.1. Symmetric pole-on disk
For a pole-on disk, $i_\text{tilt}=0^\circ$.
For a symmetric disk, $e=0$ and the position angle (pa) and argument of pericenter ($\omega$) have no impact.
We choose a semi-major axis of 70 a.u., a vertical profile with a gaussian distribution ($\gamma=2$), a reference scale height of 3 a.u. at the semi-major axis of the disk, and inner and outer exponent $\alpha_{in}=12$ and $\alpha_{out}=-12$
End of explanation
fake_disk1 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws','ain':alpha_in,'aout':alpha_out,
                                              'a':a,'e':0.0,'ksi0':ksi0,'gamma':gamma,'beta':beta},
                                spf_dico={'name':'HG', 'g':0., 'polar':False}, 
                                flux_max=1.)
Explanation: Then create your disk model
End of explanation
fake_disk1_map = fake_disk1.compute_scattered_light()
plot_frames(fake_disk1_map, grid=False, size_factor=6)
Explanation: The method compute_scattered_light returns the synthetic image of the disk.
End of explanation
fake_disk1.print_info()
Explanation: You can print some info on the geometrical properties of the model, the dust distribution parameters, the numerical integration parameters and the phase function parameters (detailed later). 
This can be useful because, in addition to reminding all the parameters used in the model, it also computes some properties such as the radial FWHM of the disk.
End of explanation
fake_disk1 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':-3,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'HG', 'g':0., 'polar':False},
                                flux_max=1.)
fake_disk1_map = fake_disk1.compute_scattered_light()
plot_frames(fake_disk1_map, grid=False, size_factor=6)
fake_disk1.print_info()
Explanation: As a side note, if $\alpha_{in} \ne \alpha_{out}$, then the peak surface density of the disk is not located at the reference radius $a$.
End of explanation
itilt = 76 # inclination of your disk in degreess
fake_disk2 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'HG', 'g':0., 'polar':False},
                                flux_max=1.)
fake_disk2_map = fake_disk2.compute_scattered_light()
plot_frames(fake_disk2_map, grid=False, size_factor=6)
Explanation: 6.2.2. Inclined symmetric disk
End of explanation
fake_disk2 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta,
                                              'dens_at_r0':1e6},
                                spf_dico={'name':'HG', 'g':0, 'polar':False})
fake_disk2_map = fake_disk2.compute_scattered_light()
plot_frames(fake_disk2_map, grid=False, size_factor=6)
Explanation: The position angle of the disk is 0 (e.g. north).  The phase function is asymmetric, the reason why the north and south ansae appear brighter is because the disk is not flat: it has a certain scale height and there is more dust intercepted along the line of sight in the ansae.
Note that we decided here to normalize the disk to a maximum brightness of 1, using the option flux_max=1.. This is not the only option available and you can decide to paramterize $\rho_0$ instead, using the keyword dens_at_r0 which directly specifies $\rho_0$.
End of explanation
fake_disk2 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=90, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':2, 'gamma':gamma, 'beta':beta, 
                                              'dens_at_r0':1e6},
                                spf_dico={'name':'HG', 'g':0, 'polar':False})
fake_disk2_map = fake_disk2.compute_scattered_light()
plot_frames(fake_disk2_map, grid=False, size_factor=6)
Explanation: Warning ! The code does not handle perfectly edge-on disks. There is a maximum inclination close to edge-on beyond which it cannot create an image. In practice this is not a limitation as the convolution by the PSF always makes it impossible to disentangle between a close to edge-on disk and a perfectly edge-on disk.
End of explanation
g=0.4
fake_disk3 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'HG', 'g':g, 'polar':False},
                                flux_max=1.)
Explanation: 6.2.3. Inclined symmetric disk with anisotropy of scattering
6.2.3.1. Simple Henyey-Greenstein phase function
We parametrize the phase function by a Henyey Greenstein phase function, with an asymmetry parameter g. An isotropic phase function has $g=0$, forward scattering is represented by $0<g\leq1$ and backward scattering is represented by $-1\leq g<0$
End of explanation
fake_disk3.phase_function.plot_phase_function()
fake_disk3_map = fake_disk3.compute_scattered_light()
plot_frames(fake_disk3_map, grid=False, size_factor=6)
Explanation: You can plot how the phase function look like:
End of explanation
g1=0.6
g2=-0.4
weight1=0.7
fake_disk4 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'DoubleHG', 'g':[g1,g2], 'weight':weight1,
                                          'polar':False},
                                flux_max=1)
fake_disk4.phase_function.plot_phase_function()
fake_disk4_map = fake_disk4.compute_scattered_light()
plot_frames(fake_disk4_map, grid=False, size_factor=6)
Explanation: The forward side is brighter. 
6.2.3.2. Double Henyey-Greenstein phase function
A double Henyey Greenstein (HG) phase function is simply a linear combination of 2 simple HG phase function. It is therefore parametrized by $g_1$ and $g_2$ the 2 asymmetry parameters of each HG and the weight (between 0 and 1) of the first HG phase function. Typically a double HG is used to represent a combination of forward scattering ($g_1>0$) and backward scattering ($g_2<1$)
End of explanation
kind='cubic' #kind must be either "linear", "nearest", "zero", "slinear", "quadratic" or "cubic"
spf_dico = dict({'phi':[0,  60,  90, 120, 180],
                 'spf':[1, 0.4, 0.3, 0.3, 0.5],
                 'name':'interpolated', 'polar':False, 'kind':kind})
fake_disk5 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico=spf_dico, flux_max=1)
fake_disk5.phase_function.plot_phase_function()
fake_disk5_map = fake_disk5.compute_scattered_light()
plot_frames(fake_disk5_map, grid=False, size_factor=6)
Explanation: 6.2.3.3. Custom phase function
In some cases, a HG phase function (simple or double) cannot represent well the behaviour of the dust. The code is modular and you can propose new prescriptions for the phase functions if you need, or you can also create a custom phase function.
End of explanation
fake_disk6 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=0, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':0.0, 'ksi0':ksi0, 'gamma':gamma, 
                                              'beta':beta, 'dens_at_r0':1e6},
                                spf_dico={'name':'HG', 'g':0, 'polar':True})
fake_disk6.phase_function.plot_phase_function()
fake_disk6_map = fake_disk6.compute_scattered_light()
plot_frames(fake_disk6_map, grid=False, size_factor=6)
Explanation: 6.2.3.4. Representing a polarised phase function
If you are trying to reproduce the polarised intensity of a disk (for instance Stokes $Q_\phi$ image), you may want to add on top of the scattering phase function, a modulation representing the degree of linear polarisation.
This can be done by setting the polar keyword to True and in this case, the model assumes a Rayleigh-like degree of linear polarisation parametrized by $(1-(\cos \phi)^2) / (1+(\cos \phi)^2)$ where $\phi$ is the scattering angle.
End of explanation
e=0.4 # eccentricity in degrees
omega=30 # argument of pericenter
fake_disk7 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=0, omega=omega, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':e, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'HG', 'g':g, 'polar':False},
                                flux_max=1.)
fake_disk7_map = fake_disk7.compute_scattered_light()
plot_frames(fake_disk7_map, grid=False, size_factor=6)
Explanation: You can combine this Rayleigh-like degree of linear polarisation with any phase function (simple HG, double HG or custom type).
6.2.4. Asymmetric disk
Be careful here !
There is no consensus in the community on how to parametrize an eccentric dust distribution, so keep in mind that the convention described in section 1.2 is only one way to do so, but does not mean the dust density distribution in an eccentric disk follows this prescription. For instance, around the pericenter particle velocities are higher and one expects more collision to happen which can create an overdensity of particles compared to other regions of the disk. Conversely, particles stay longer at the apocenter because of Kepler's third law, which means that one could also expect a higher density at apocenter... All these physical phenomena are not described in this model.
Let's start woth a pole-on disk to be insensitive to phase function effects
End of explanation
fake_disk7 = ScatteredLightDisk(nx=nx, ny=ny, distance=dstar,
                                itilt=itilt, omega=omega, pxInArcsec=pixel_scale, pa=0,
                                density_dico={'name':'2PowerLaws', 'ain':alpha_in, 'aout':alpha_out,
                                              'a':a, 'e':e, 'ksi0':ksi0, 'gamma':gamma, 'beta':beta},
                                spf_dico={'name':'HG', 'g':g, 'polar':False},
                                flux_max=1.)
fake_disk7_map = fake_disk7.compute_scattered_light()
plot_frames(fake_disk7_map, grid=False, size_factor=6)
Explanation: The brightness asymmetry here is entirely due to the fact that the brightness at one point in the disk is inversely proportional to the squared distance to the star.
Once you incline the disk, you start seeing the competing effect of the phase function and eccentricity.
End of explanation
plot_frames(fake_disk3_map, grid=False, size_factor=6)
nframes = 30
# we assume we have 60º of parallactic angle rotation centered around meridian
parang_amplitude = 60
derotation_angles = np.linspace(-parang_amplitude/2, parang_amplitude/2, nframes)
start = time_ini()
cube_fake_disk3 = cube_inject_fakedisk(fake_disk3_map, -derotation_angles, imlib='vip-fft')
timing(start)
Explanation: 6.3. Forward modeling of disks
Let's start from our inclined simple HG symmeric disk fake_disk3_map and assume we observe this disk as part of an ADI sequence of 30 images
End of explanation
cube_fake_disk3.shape
Explanation: cube_fake_disk3 is now a cube of 30 frames, where the disk has been injected at the correct position angle.
End of explanation
plot_frames((cube_fake_disk3[0], cube_fake_disk3[nframes//2], cube_fake_disk3[nframes-1]), 
            grid=False, size_factor=3)
Explanation: Let's visualize the first, middle and last image of the cube.
End of explanation
cadi_fake_disk3 = median_sub(cube_fake_disk3, derotation_angles, imlib='vip-fft')
plot_frames((fake_disk3_map, cadi_fake_disk3), grid=False, size_factor=4)
Explanation: We can now process this cube with median-ADI for instance:
End of explanation
psf = create_synth_psf(model='gauss', shape=(11, 11), fwhm=4.)
plot_frames(psf, grid=True, size_factor=2)
Explanation: The example above shows a typical bias that can be induced by ADI on extended disk signals (Milli et al. 2012).
So far we have not dealt with convolution effects. In practice the image of a disk is convolved by the instrumental PSF. 
Let's assume here an instrument having a gaussian PSF with FWHM = 4px, and create a synthetic PSF using the create_synth_psf function:
End of explanation
cube_fake_disk3_convolved = cube_inject_fakedisk(fake_disk3_map, -derotation_angles,
                                                 psf=psf, imlib='vip-fft')
cadi_fake_disk3_convolved = median_sub(cube_fake_disk3_convolved, derotation_angles, imlib='vip-fft')
plot_frames((fake_disk3_map, cadi_fake_disk3, cadi_fake_disk3_convolved), grid=False, size_factor=4)
Explanation: Then we inject the disk in the cube and convolve each frame by the PSF
End of explanation
<END_TASK>