diff --git a/Makefile b/Makefile index c1915502..c489b22f 100644 --- a/Makefile +++ b/Makefile @@ -18,8 +18,8 @@ pypi-test: docs: rm -rf docs/api - sphinx-apidoc -o docs vip_hci - cd docs/ + sphinx-apidoc -o docs/source vip_hci + cd docs/source/ $(MAKE) -C docs clean $(MAKE) -C docs html diff --git a/docs/source/vip_hci.andromeda.rst b/docs/source/vip_hci.andromeda.rst deleted file mode 100644 index 0e4af1c5..00000000 --- a/docs/source/vip_hci.andromeda.rst +++ /dev/null @@ -1,44 +0,0 @@ -vip\_hci.andromeda package -========================== - -.. automodule:: vip_hci.andromeda - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.andromeda.andromeda module ------------------------------------ - -.. automodule:: vip_hci.andromeda.andromeda - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.andromeda.fit module ------------------------------ - -.. automodule:: vip_hci.andromeda.fit - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.andromeda.shift module -------------------------------- - -.. automodule:: vip_hci.andromeda.shift - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.andromeda.utils module -------------------------------- - -.. automodule:: vip_hci.andromeda.utils - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.conf.rst b/docs/source/vip_hci.conf.rst deleted file mode 100644 index 5c86ae79..00000000 --- a/docs/source/vip_hci.conf.rst +++ /dev/null @@ -1,44 +0,0 @@ -vip\_hci.conf package -===================== - -.. automodule:: vip_hci.conf - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.conf.mem module ------------------------- - -.. automodule:: vip_hci.conf.mem - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.conf.param module --------------------------- - -.. automodule:: vip_hci.conf.param - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.conf.timing module ---------------------------- - -.. automodule:: vip_hci.conf.timing - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.conf.utils\_conf module --------------------------------- - -.. automodule:: vip_hci.conf.utils_conf - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.config.rst b/docs/source/vip_hci.config.rst new file mode 100644 index 00000000..54a35409 --- /dev/null +++ b/docs/source/vip_hci.config.rst @@ -0,0 +1,45 @@ +vip\_hci.config package +======================= + +Submodules +---------- + +vip\_hci.config.mem module +-------------------------- + +.. automodule:: vip_hci.config.mem + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.config.param module +---------------------------- + +.. automodule:: vip_hci.config.param + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.config.timing module +----------------------------- + +.. automodule:: vip_hci.config.timing + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.config.utils\_conf module +---------------------------------- + +.. automodule:: vip_hci.config.utils_conf + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: vip_hci.config + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.exlib.rst b/docs/source/vip_hci.exlib.rst deleted file mode 100644 index 681d8173..00000000 --- a/docs/source/vip_hci.exlib.rst +++ /dev/null @@ -1,20 +0,0 @@ -vip\_hci.exlib package -====================== - -.. automodule:: vip_hci.exlib - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.exlib.iuwt module --------------------------- - -.. automodule:: vip_hci.exlib.iuwt - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.fits.rst b/docs/source/vip_hci.fits.rst index 87ce21b7..350c68fd 100644 --- a/docs/source/vip_hci.fits.rst +++ b/docs/source/vip_hci.fits.rst @@ -1,11 +1,6 @@ vip\_hci.fits package ===================== -.. automodule:: vip_hci.fits - :members: - :undoc-members: - :show-inheritance: - Submodules ---------- @@ -13,8 +8,14 @@ vip\_hci.fits.fits module ------------------------- .. automodule:: vip_hci.fits.fits - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: +Module contents +--------------- +.. automodule:: vip_hci.fits + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.fm.rst b/docs/source/vip_hci.fm.rst new file mode 100644 index 00000000..fd4a4ec4 --- /dev/null +++ b/docs/source/vip_hci.fm.rst @@ -0,0 +1,93 @@ +vip\_hci.fm package +=================== + +Submodules +---------- + +vip\_hci.fm.fakecomp module +--------------------------- + +.. automodule:: vip_hci.fm.fakecomp + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.fakedisk module +--------------------------- + +.. automodule:: vip_hci.fm.fakedisk + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.negfc\_fmerit module +-------------------------------- + +.. automodule:: vip_hci.fm.negfc_fmerit + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.negfc\_mcmc module +------------------------------ + +.. automodule:: vip_hci.fm.negfc_mcmc + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.negfc\_nested module +-------------------------------- + +.. automodule:: vip_hci.fm.negfc_nested + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.negfc\_simplex module +--------------------------------- + +.. automodule:: vip_hci.fm.negfc_simplex + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.negfc\_speckle\_noise module +---------------------------------------- + +.. automodule:: vip_hci.fm.negfc_speckle_noise + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.scattered\_light\_disk module +----------------------------------------- + +.. automodule:: vip_hci.fm.scattered_light_disk + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.utils\_mcmc module +------------------------------ + +.. automodule:: vip_hci.fm.utils_mcmc + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.fm.utils\_negfc module +------------------------------- + +.. automodule:: vip_hci.fm.utils_negfc + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: vip_hci.fm + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.frdiff.rst b/docs/source/vip_hci.frdiff.rst deleted file mode 100644 index 5b982b88..00000000 --- a/docs/source/vip_hci.frdiff.rst +++ /dev/null @@ -1,20 +0,0 @@ -vip\_hci.frdiff package -======================= - -.. automodule:: vip_hci.frdiff - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.frdiff.framediff module --------------------------------- - -.. automodule:: vip_hci.frdiff.framediff - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.invprob.rst b/docs/source/vip_hci.invprob.rst new file mode 100644 index 00000000..34f60ba5 --- /dev/null +++ b/docs/source/vip_hci.invprob.rst @@ -0,0 +1,29 @@ +vip\_hci.invprob package +======================== + +Submodules +---------- + +vip\_hci.invprob.andromeda module +--------------------------------- + +.. automodule:: vip_hci.invprob.andromeda + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.invprob.utils\_andro module +------------------------------------ + +.. automodule:: vip_hci.invprob.utils_andro + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: vip_hci.invprob + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.leastsq.rst b/docs/source/vip_hci.leastsq.rst deleted file mode 100644 index 4e3619c2..00000000 --- a/docs/source/vip_hci.leastsq.rst +++ /dev/null @@ -1,20 +0,0 @@ -vip\_hci.leastsq package -======================== - -.. automodule:: vip_hci.leastsq - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.leastsq.leastsq module -------------------------------- - -.. automodule:: vip_hci.leastsq.leastsq - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.llsg.rst b/docs/source/vip_hci.llsg.rst deleted file mode 100644 index 2586593f..00000000 --- a/docs/source/vip_hci.llsg.rst +++ /dev/null @@ -1,28 +0,0 @@ -vip\_hci.llsg package -===================== - -.. automodule:: vip_hci.llsg - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.llsg.llsg module -------------------------- - -.. automodule:: vip_hci.llsg.llsg - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.llsg.thresholding module ---------------------------------- - -.. automodule:: vip_hci.llsg.thresholding - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.medsub.rst b/docs/source/vip_hci.medsub.rst deleted file mode 100644 index 69561507..00000000 --- a/docs/source/vip_hci.medsub.rst +++ /dev/null @@ -1,20 +0,0 @@ -vip\_hci.medsub package -======================= - -.. automodule:: vip_hci.medsub - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.medsub.medsub\_source module -------------------------------------- - -.. automodule:: vip_hci.medsub.medsub_source - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.metrics.rst b/docs/source/vip_hci.metrics.rst index 9e729eda..6864e591 100644 --- a/docs/source/vip_hci.metrics.rst +++ b/docs/source/vip_hci.metrics.rst @@ -1,11 +1,6 @@ vip\_hci.metrics package ======================== -.. automodule:: vip_hci.metrics - :members: - :undoc-members: - :show-inheritance: - Submodules ---------- @@ -13,80 +8,46 @@ vip\_hci.metrics.contrcurve module ---------------------------------- .. automodule:: vip_hci.metrics.contrcurve - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.metrics.detection module --------------------------------- .. automodule:: vip_hci.metrics.detection - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.metrics.dust\_distribution module ------------------------------------------- - -.. automodule:: vip_hci.metrics.dust_distribution - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.metrics.fakecomp module --------------------------------- - -.. automodule:: vip_hci.metrics.fakecomp - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.metrics.fakedisk module --------------------------------- - -.. automodule:: vip_hci.metrics.fakedisk - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.metrics.phase\_function module ---------------------------------------- - -.. automodule:: vip_hci.metrics.phase_function - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.metrics.roc module --------------------------- .. automodule:: vip_hci.metrics.roc - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.metrics.scattered\_light\_disk module ----------------------------------------------- - -.. automodule:: vip_hci.metrics.scattered_light_disk - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.metrics.snr\_source module ----------------------------------- .. automodule:: vip_hci.metrics.snr_source - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.metrics.stim module ---------------------------- .. automodule:: vip_hci.metrics.stim - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: +Module contents +--------------- +.. automodule:: vip_hci.metrics + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.negfc.rst b/docs/source/vip_hci.negfc.rst deleted file mode 100644 index 295963e5..00000000 --- a/docs/source/vip_hci.negfc.rst +++ /dev/null @@ -1,52 +0,0 @@ -vip\_hci.negfc package -====================== - -.. automodule:: vip_hci.negfc - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.negfc.mcmc\_sampling module ------------------------------------- - -.. automodule:: vip_hci.negfc.mcmc_sampling - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.negfc.nested\_sampling module --------------------------------------- - -.. automodule:: vip_hci.negfc.nested_sampling - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.negfc.simplex\_fmerit module -------------------------------------- - -.. automodule:: vip_hci.negfc.simplex_fmerit - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.negfc.simplex\_optim module ------------------------------------- - -.. automodule:: vip_hci.negfc.simplex_optim - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.negfc.utils\_negfc module ----------------------------------- - -.. automodule:: vip_hci.negfc.utils_negfc - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.nmf.rst b/docs/source/vip_hci.nmf.rst deleted file mode 100644 index 0af3aa88..00000000 --- a/docs/source/vip_hci.nmf.rst +++ /dev/null @@ -1,20 +0,0 @@ -vip\_hci.nmf package -==================== - -.. automodule:: vip_hci.nmf - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.nmf.nmf\_fullfr module -------------------------------- - -.. automodule:: vip_hci.nmf.nmf_fullfr - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.pca.rst b/docs/source/vip_hci.pca.rst deleted file mode 100644 index e8443464..00000000 --- a/docs/source/vip_hci.pca.rst +++ /dev/null @@ -1,44 +0,0 @@ -vip\_hci.pca package -==================== - -.. automodule:: vip_hci.pca - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.pca.pca\_fullfr module -------------------------------- - -.. automodule:: vip_hci.pca.pca_fullfr - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.pca.pca\_local module ------------------------------- - -.. automodule:: vip_hci.pca.pca_local - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.pca.svd module ------------------------ - -.. automodule:: vip_hci.pca.svd - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.pca.utils\_pca module ------------------------------- - -.. automodule:: vip_hci.pca.utils_pca - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.preproc.rst b/docs/source/vip_hci.preproc.rst index e57aca6f..991b96c0 100644 --- a/docs/source/vip_hci.preproc.rst +++ b/docs/source/vip_hci.preproc.rst @@ -1,11 +1,6 @@ vip\_hci.preproc package ======================== -.. automodule:: vip_hci.preproc - :members: - :undoc-members: - :show-inheritance: - Submodules ---------- @@ -13,72 +8,78 @@ vip\_hci.preproc.badframes module --------------------------------- .. automodule:: vip_hci.preproc.badframes - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.badpixremoval module ------------------------------------- .. automodule:: vip_hci.preproc.badpixremoval - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.cosmetics module --------------------------------- .. automodule:: vip_hci.preproc.cosmetics - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.derotation module ---------------------------------- .. automodule:: vip_hci.preproc.derotation - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.parangles module --------------------------------- .. automodule:: vip_hci.preproc.parangles - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.recentering module ----------------------------------- .. automodule:: vip_hci.preproc.recentering - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.rescaling module --------------------------------- .. automodule:: vip_hci.preproc.rescaling - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.skysubtraction module -------------------------------------- .. automodule:: vip_hci.preproc.skysubtraction - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.preproc.subsampling module ----------------------------------- .. automodule:: vip_hci.preproc.subsampling - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: +Module contents +--------------- +.. automodule:: vip_hci.preproc + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.psfsub.rst b/docs/source/vip_hci.psfsub.rst new file mode 100644 index 00000000..330caacd --- /dev/null +++ b/docs/source/vip_hci.psfsub.rst @@ -0,0 +1,93 @@ +vip\_hci.psfsub package +======================= + +Submodules +---------- + +vip\_hci.psfsub.framediff module +-------------------------------- + +.. automodule:: vip_hci.psfsub.framediff + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.llsg module +--------------------------- + +.. automodule:: vip_hci.psfsub.llsg + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.loci module +--------------------------- + +.. automodule:: vip_hci.psfsub.loci + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.medsub module +----------------------------- + +.. automodule:: vip_hci.psfsub.medsub + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.nmf\_fullfr module +---------------------------------- + +.. automodule:: vip_hci.psfsub.nmf_fullfr + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.nmf\_local module +--------------------------------- + +.. automodule:: vip_hci.psfsub.nmf_local + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.pca\_fullfr module +---------------------------------- + +.. automodule:: vip_hci.psfsub.pca_fullfr + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.pca\_local module +--------------------------------- + +.. automodule:: vip_hci.psfsub.pca_local + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.svd module +-------------------------- + +.. automodule:: vip_hci.psfsub.svd + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.psfsub.utils\_pca module +--------------------------------- + +.. automodule:: vip_hci.psfsub.utils_pca + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: vip_hci.psfsub + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.rst b/docs/source/vip_hci.rst index bb4c9c4e..86d2ad1a 100644 --- a/docs/source/vip_hci.rst +++ b/docs/source/vip_hci.rst @@ -11,20 +11,13 @@ Subpackages .. toctree:: - vip_hci.andromeda - vip_hci.conf - vip_hci.exlib + vip_hci.config vip_hci.fits - vip_hci.frdiff - vip_hci.leastsq - vip_hci.llsg - vip_hci.medsub + vip_hci.fm + vip_hci.invprob vip_hci.metrics - vip_hci.negfc - vip_hci.nmf - vip_hci.pca vip_hci.preproc - vip_hci.specfit + vip_hci.psfsub vip_hci.stats vip_hci.var diff --git a/docs/source/vip_hci.specfit.rst b/docs/source/vip_hci.specfit.rst deleted file mode 100644 index 3fa4dd31..00000000 --- a/docs/source/vip_hci.specfit.rst +++ /dev/null @@ -1,52 +0,0 @@ -vip\_hci.specfit package -======================== - -.. automodule:: vip_hci.specfit - :members: - :undoc-members: - :show-inheritance: - -Submodules ----------- - -vip\_hci.specfit.chi module ---------------------------- - -.. automodule:: vip_hci.specfit.chi - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.specfit.mcmc\_sampling\_spec module --------------------------------------------- - -.. automodule:: vip_hci.specfit.mcmc_sampling_spec - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.specfit.model\_resampling module ------------------------------------------ - -.. automodule:: vip_hci.specfit.model_resampling - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.specfit.template\_lib module -------------------------------------- - -.. automodule:: vip_hci.specfit.template_lib - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.specfit.utils\_spec module ------------------------------------ - -.. automodule:: vip_hci.specfit.utils_spec - :members: - :undoc-members: - :show-inheritance: - - diff --git a/docs/source/vip_hci.stats.rst b/docs/source/vip_hci.stats.rst index 2e724c49..ad5318a8 100644 --- a/docs/source/vip_hci.stats.rst +++ b/docs/source/vip_hci.stats.rst @@ -1,44 +1,53 @@ vip\_hci.stats package ====================== -.. automodule:: vip_hci.stats - :members: - :undoc-members: - :show-inheritance: - Submodules ---------- +vip\_hci.stats.bkg\_proba module +-------------------------------- + +.. automodule:: vip_hci.stats.bkg_proba + :members: + :undoc-members: + :show-inheritance: + vip\_hci.stats.clip\_sigma module --------------------------------- .. automodule:: vip_hci.stats.clip_sigma - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.stats.distances module ------------------------------- .. automodule:: vip_hci.stats.distances - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.stats.im\_stats module ------------------------------- .. automodule:: vip_hci.stats.im_stats - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.stats.utils\_stats module ---------------------------------- .. automodule:: vip_hci.stats.utils_stats - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: +Module contents +--------------- +.. automodule:: vip_hci.stats + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/vip_hci.var.rst b/docs/source/vip_hci.var.rst index 3c42f717..fc7fc905 100644 --- a/docs/source/vip_hci.var.rst +++ b/docs/source/vip_hci.var.rst @@ -1,11 +1,6 @@ vip\_hci.var package ==================== -.. automodule:: vip_hci.var - :members: - :undoc-members: - :show-inheritance: - Submodules ---------- @@ -13,38 +8,46 @@ vip\_hci.var.coords module -------------------------- .. automodule:: vip_hci.var.coords - :members: - :undoc-members: - :show-inheritance: - + :members: + :undoc-members: + :show-inheritance: + vip\_hci.var.filters module --------------------------- .. automodule:: vip_hci.var.filters - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: vip\_hci.var.fit\_2d module --------------------------- .. automodule:: vip_hci.var.fit_2d - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: + +vip\_hci.var.iuwt module +------------------------ + +.. automodule:: vip_hci.var.iuwt + :members: + :undoc-members: + :show-inheritance: vip\_hci.var.shapes module -------------------------- .. automodule:: vip_hci.var.shapes - :members: - :undoc-members: - :show-inheritance: - -vip\_hci.var.iuwt module ---------------------------- + :members: + :undoc-members: + :show-inheritance: -.. automodule:: vip_hci.var.iuwt - :members: - :undoc-members: - :show-inheritance: \ No newline at end of file +Module contents +--------------- + +.. automodule:: vip_hci.var + :members: + :undoc-members: + :show-inheritance: diff --git a/readthedocs.yml b/readthedocs.yml index f2f64853..2d1e821a 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -1,5 +1,8 @@ # .readthedocs.yml +sphinx: + configuration: docs/source/conf.py + build: image: latest diff --git a/setup.py b/setup.py index 74b59836..7eb45b53 100644 --- a/setup.py +++ b/setup.py @@ -60,9 +60,9 @@ def resource(*args): PACKAGES = ['vip_hci', 'vip_hci.config', 'vip_hci.fits', + 'vip_hci.fm', 'vip_hci.invprob', 'vip_hci.metrics', - 'vip_hci.negfc', 'vip_hci.preproc', 'vip_hci.psfsub', 'vip_hci.stats', diff --git a/tests/test_metrics_fakecomp.py b/tests/test_fm_fakecomp.py similarity index 96% rename from tests/test_metrics_fakecomp.py rename to tests/test_fm_fakecomp.py index cf04032d..b46daad9 100644 --- a/tests/test_metrics_fakecomp.py +++ b/tests/test_fm_fakecomp.py @@ -1,10 +1,10 @@ """ -Tests for metrics/fakecomp.py +Tests for fm/fakecomp.py """ from .helpers import aarc, np, param, parametrize, fixture, filterwarnings -from vip_hci.metrics.fakecomp import cube_inject_companions, normalize_psf +from vip_hci.fm import cube_inject_companions, normalize_psf # ===== utility functions diff --git a/tests/test_fm_fakedisk.py b/tests/test_fm_fakedisk.py new file mode 100644 index 00000000..f03e22fa --- /dev/null +++ b/tests/test_fm_fakedisk.py @@ -0,0 +1,93 @@ +""" +Tests for fm/fakecomp.py + +""" + +from .helpers import aarc, np, fixture +from vip_hci.fm import cube_inject_fakedisk, cube_inject_trace + + +# ===== utility functions + +@fixture(scope="module", params=["3D"]) +def dataset(request): + """ + Create 3D and 4D datasets for use with ``test_cube_inject_companions``. + + """ + if request.param == "3D": + cube = np.zeros((3,25,25)) + psf = np.ones((1,1)) + + angles = np.array([0, 90, 180]) + + return cube, psf, angles + + +def test_cube_inject_fakedisk(dataset): + """ + Verify position of injected disk image with 1 value, for 3D and 4D cases. + """ + def _expected(): + """ + Expected positions. + """ + return [(15, 12), (12, 9), (9, 12)] + + psf = np.zeros((25, 25)) + psf[15,12] = 1 + + _, _, angles = dataset + + cube = cube_inject_fakedisk(psf, angle_list=angles) + + # find coords + coords = [] + for i in range(cube.shape[0]): + max_idx = np.argmax(cube[i]) + coords.append(np.unravel_index(max_idx, cube[0].shape)) + + yx_expected = _expected() + + aarc(coords, yx_expected) + + +def test_cube_inject_trace(dataset): + """ + Verify position of injected disk image with 1 value, for 3D and 4D cases. + """ + def _expected(ang): + """ + Expected positions. + """ + if ang == 0: + return [(7, 12), (12, 8), (15, 12)] + elif ang == 90: + return [(12, 7), (12, 15), (16, 12)] + elif ang == 180: + return [(9, 12), (12, 16), (17, 12)] + + cube, psf, angles = dataset + + rads = [3,4,5] + thetas = [90,180,270] + + cube = cube_inject_trace(cube, psf, angles, flevel=1, + rad_dists=rads, theta=thetas, + plsc=0.01225, n_branches=1, imlib='vip-fft', + interpolation='lanczos4', verbose=True) + + for i in range(cube.shape[0]): + # find coords of trace in each image of the cube + coords = [] + nspi = len(rads) + frame_tmp = cube[i].copy() + for s in range(nspi): + max_idx = np.argmax(frame_tmp) + coords_tmp = np.unravel_index(max_idx, frame_tmp.shape) + coords.append(coords_tmp) + frame_tmp[coords_tmp]=0 + + yx_expected = _expected(angles[i]) + + aarc(coords, yx_expected) diff --git a/tests/test_fm_scatteredlightdisk.py b/tests/test_fm_scatteredlightdisk.py new file mode 100644 index 00000000..c2d00738 --- /dev/null +++ b/tests/test_fm_scatteredlightdisk.py @@ -0,0 +1,36 @@ +""" +Tests for fm/scattered_light_disk.py + +""" + +from .helpers import aarc, parametrize, param +from vip_hci.fm.scattered_light_disk import DustEllipticalDistribution2PowerLaws + + +@parametrize("r", + [ + param(55), + param(60), + param(65) + ]) +def test_dust_distribution(r): + """ + Verify dust density at different radii, for a distribution centered at a=60 + """ + def t_expected(r): + """ + Expected positions. + """ + if r == 55: + return 0.844237 + elif r == 60: + return 1 + elif r == 65: + return 0.864574 + test = DustEllipticalDistribution2PowerLaws() + test.set_radial_density(ain=5., aout=-5., a=60., e=0.) + costheta = 0. + z = 0. + t = test.density_cylindrical(r, costheta, z) + + aarc(t, t_expected(r)) \ No newline at end of file diff --git a/tests/test_metrics_snr.py b/tests/test_metrics_snr.py index e25e4c86..ac66e4c0 100644 --- a/tests/test_metrics_snr.py +++ b/tests/test_metrics_snr.py @@ -65,7 +65,7 @@ def test_snrmap_fast(get_frame): y1, x1 = np.where(snmap == snmap.max()) assert np.allclose(x1, x0, atol=atol) and np.allclose(y1, y0, atol=atol) - + def test_frame_report(get_frame): frame, positions = get_frame y0, x0 = positions diff --git a/tests/test_metrics_stim.py b/tests/test_metrics_stim.py new file mode 100644 index 00000000..420733bb --- /dev/null +++ b/tests/test_metrics_stim.py @@ -0,0 +1,59 @@ +""" +Tests for metrics/stim.py + +""" + +import copy +from .helpers import fixture, np +from vip_hci.psfsub import pca +from vip_hci.hci_dataset import Frame +from vip_hci.metrics import stim_map, inverse_stim_map + + +@fixture(scope="module") +def get_frame(example_dataset_adi): + """ + Inject a fake companion into an example cube. + + Parameters + ---------- + example_dataset_adi : fixture + Taken automatically from ``conftest.py``. + + Returns + ------- + frame : VIP Frame + planet_position : tuple(y, x) + + """ + dsi = copy.copy(example_dataset_adi) + # we chose a shallow copy, as we will not use any in-place operations + # (like +=). Using `deepcopy` would be safer, but consume more memory. + + print("producing a final frame...") + res = pca(dsi.cube, dsi.angles, ncomp=10, full_output=True) + res_frame = res[0] + res_cube = res[-2] + res_der_cube = res[-1] + frame = Frame(res_frame, fwhm=dsi.fwhm) + return frame, (63, 63), res_cube, res_der_cube, dsi.angles + + +atol = 2 +plot = False + +def test_stimmap(get_frame): + frame, positions, _, res_der_cube, _ = get_frame + y0, x0 = positions + stimap = stim_map(res_der_cube) + y1, x1 = np.where(stimap == stimap.max()) + assert np.allclose(x1, x0, atol=atol) and np.allclose(y1, y0, atol=atol) + +def test_normstimmap(get_frame): + frame, positions, res_cube, res_der_cube, angles = get_frame + y0, x0 = positions + stimap = stim_map(res_der_cube) + inv_stimap = inverse_stim_map(res_cube, angles) + norm_stimap = stimap/np.amax(inv_stimap) + y1, x1 = np.where(norm_stimap == norm_stimap.max()) + assert np.allclose(x1, x0, atol=atol) and np.allclose(y1, y0, atol=atol) diff --git a/vip_hci/__init__.py b/vip_hci/__init__.py index 288e8dc4..d78fd15a 100644 --- a/vip_hci/__init__.py +++ b/vip_hci/__init__.py @@ -5,7 +5,7 @@ from . import fits from . import invprob from . import psfsub -from . import negfc +from . import fm from . import metrics from . import stats from . import var diff --git a/vip_hci/fm/__init__.py b/vip_hci/fm/__init__.py new file mode 100644 index 00000000..b8fc57ef --- /dev/null +++ b/vip_hci/fm/__init__.py @@ -0,0 +1,28 @@ +r""" +Subpackge ``fm`` contains an ensemble of algorithms for forward modeling, +including scattered light disk model creation, fake planet and fake disk +injection, and planet position+flux estimation through the negative fake +companion algorithm (NEGFC). NEGFC can work with either a simplex (Nelder-Mead) +minimizer or coupled with Monte Carlo methods for posterior sampling. The +latter allows the estimation of uncertainties for the photometry and position +of the companion. Two possible ways of sampling the posteriors are implemented: +using ``emcee`` and its Affine Invariant MCMC or ``nestle`` with either a +single or multiple ellipsoid nested sampling procedure. + +The main idea of the NegFC is to inject negative fake companions (candidates) +with varying position and flux in the original cube and minimize a figure of +merit corresponding to the residuals in an aperture at the initial estimate for +the location of the companion, in the post-processed image. + +""" + +from .scattered_light_disk import * +from .fakecomp import * +from .fakedisk import * +from .negfc_fmerit import * +from .negfc_mcmc import * +from .negfc_nested import * +from .negfc_simplex import * +from .negfc_speckle_noise import * +from .utils_mcmc import * +from .utils_negfc import * \ No newline at end of file diff --git a/vip_hci/metrics/fakecomp.py b/vip_hci/fm/fakecomp.py similarity index 98% rename from vip_hci/metrics/fakecomp.py rename to vip_hci/fm/fakecomp.py index fd5faf78..a4e08457 100644 --- a/vip_hci/metrics/fakecomp.py +++ b/vip_hci/fm/fakecomp.py @@ -14,7 +14,14 @@ import numpy as np from scipy import stats from scipy.interpolate import interp1d +from packaging import version import photutils +if version.parse(photutils.__version__) >= version.parse('0.3'): + # for photutils version >= '0.3' use photutils.centroids.centroid_com + from photutils.centroids import centroid_com as cen_com +else: + # for photutils version < '0.3' use photutils.centroid_com + import photutils.centroid_com as cen_com from ..preproc import (cube_crop_frames, frame_shift, frame_crop, cube_shift, frame_rotate) from ..var import (frame_center, fit_2dgaussian, fit_2dairydisk, fit_2dmoffat, @@ -149,7 +156,7 @@ def cube_inject_companions(array, psf_template, angle_list, flevel, plsc, if size_fc%2: # new convention w -= 1 sty = int(ceny) - w - sty = int(cenx) - w + stx = int(cenx) - w # fake companion in the center of a zeros frame fc_fr = np.zeros_like(array) @@ -165,7 +172,7 @@ def cube_inject_companions(array, psf_template, angle_list, flevel, plsc, else: try: for fr in range(nframes): - fc_fr[fr,sty:sty+size_fc, sty:sty+size_fc] = psf_template[fr] + fc_fr[fr,sty:sty+size_fc, stx:stx+size_fc] = psf_template[fr] except ValueError as e: print("cannot place PSF on frames. Please verify the shapes!" "psf shape: {}, array shape: {}".format(psf_template.shape, @@ -535,7 +542,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose): """ 2d case """ # we check if the psf is centered and fix it if needed cy, cx = frame_center(psf, verbose=False) - xcom, ycom = photutils.centroid_com(psf) + xcom, ycom = cen_com(psf) if not (np.allclose(cy, ycom, atol=1e-2) or np.allclose(cx, xcom, atol=1e-2)): # first we find the centroid and put it in the center of the array diff --git a/vip_hci/metrics/fakedisk.py b/vip_hci/fm/fakedisk.py similarity index 100% rename from vip_hci/metrics/fakedisk.py rename to vip_hci/fm/fakedisk.py diff --git a/vip_hci/negfc/simplex_fmerit.py b/vip_hci/fm/negfc_fmerit.py similarity index 99% rename from vip_hci/negfc/simplex_fmerit.py rename to vip_hci/fm/negfc_fmerit.py index 6084a932..ebe93155 100644 --- a/vip_hci/negfc/simplex_fmerit.py +++ b/vip_hci/fm/negfc_fmerit.py @@ -10,7 +10,7 @@ import numpy as np from hciplot import plot_frames from skimage.draw import disk -from ..metrics import cube_inject_companions +from ..fm import cube_inject_companions from ..var import (frame_center, get_annular_wedge, cube_filter_highpass, dist, get_circle) from ..psfsub import pca_annulus, pca_annular, pca diff --git a/vip_hci/negfc/mcmc_sampling.py b/vip_hci/fm/negfc_mcmc.py similarity index 99% rename from vip_hci/negfc/mcmc_sampling.py rename to vip_hci/fm/negfc_mcmc.py index 7b851bda..c25381c3 100644 --- a/vip_hci/negfc/mcmc_sampling.py +++ b/vip_hci/fm/negfc_mcmc.py @@ -22,11 +22,11 @@ from matplotlib.ticker import MaxNLocator import pickle from scipy.stats import norm -from ..metrics import cube_inject_companions +from ..fm import cube_inject_companions from ..config import time_ini, timing from ..config.utils_conf import sep from ..psfsub import pca_annulus -from .simplex_fmerit import get_values_optimize, get_mu_and_sigma +from .negfc_fmerit import get_values_optimize, get_mu_and_sigma from .utils_mcmc import gelman_rubin, autocorr_test import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) diff --git a/vip_hci/negfc/nested_sampling.py b/vip_hci/fm/negfc_nested.py similarity index 99% rename from vip_hci/negfc/nested_sampling.py rename to vip_hci/fm/negfc_nested.py index 621c71e1..7ae8b377 100644 --- a/vip_hci/negfc/nested_sampling.py +++ b/vip_hci/fm/negfc_nested.py @@ -16,7 +16,7 @@ import numpy as np from matplotlib import pyplot as plt from ..config import time_ini, timing -from .mcmc_sampling import lnlike, confidence, show_walk_plot +from .negfc_mcmc import lnlike, confidence, show_walk_plot from ..psfsub import pca_annulus def nested_negfc_sampling(init, cube, angs, plsc, psf, fwhm, annulus_width=8, diff --git a/vip_hci/negfc/simplex_optim.py b/vip_hci/fm/negfc_simplex.py similarity index 99% rename from vip_hci/negfc/simplex_optim.py rename to vip_hci/fm/negfc_simplex.py index 695225f3..6b78255e 100644 --- a/vip_hci/negfc/simplex_optim.py +++ b/vip_hci/fm/negfc_simplex.py @@ -9,7 +9,7 @@ import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize -from .simplex_fmerit import chisquare, get_mu_and_sigma +from .negfc_fmerit import chisquare, get_mu_and_sigma from ..psfsub import pca_annulus from ..var import frame_center from ..config import time_ini, timing diff --git a/vip_hci/negfc/speckle_noise.py b/vip_hci/fm/negfc_speckle_noise.py similarity index 98% rename from vip_hci/negfc/speckle_noise.py rename to vip_hci/fm/negfc_speckle_noise.py index dc17b480..1ca6df0d 100644 --- a/vip_hci/negfc/speckle_noise.py +++ b/vip_hci/fm/negfc_speckle_noise.py @@ -14,11 +14,11 @@ import matplotlib.pyplot as plt from ..config.utils_conf import pool_map, iterable #eval_func_tuple -from ..metrics import cube_inject_companions -from .simplex_optim import firstguess_simplex -from .simplex_fmerit import get_mu_and_sigma +from ..fm import cube_inject_companions +from .negfc_simplex import firstguess_simplex +from .negfc_fmerit import get_mu_and_sigma from .utils_negfc import cube_planet_free -from .mcmc_sampling import confidence +from .negfc_mcmc import confidence def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo, diff --git a/vip_hci/fm/scattered_light_disk.py b/vip_hci/fm/scattered_light_disk.py new file mode 100644 index 00000000..26b4e977 --- /dev/null +++ b/vip_hci/fm/scattered_light_disk.py @@ -0,0 +1,1112 @@ +#! /usr/bin/env python +""" +Class definition for ScatteredLightDisk, Dust_distribution and Phase_function +""" + +__author__ = 'Julien Milli' +__all__ = ['ScatteredLightDisk', + 'Dust_distribution', + 'Phase_function'] + +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import newton +from scipy.interpolate import interp1d +from ..var import frame_center + + +class ScatteredLightDisk(object): + """ + Class used to generate a synthetic disc, inspired from a light version of + the GRATER tool (GRenoble RAdiative TransfER, Augereau et al. 1999) + written originally in IDL and converted to Python by J. Milli. + """ + + def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., + pxInArcsec=0.01225, pa=0., flux_max=None, + density_dico={'name': '2PowerLaws', 'ain':5, 'aout':-5, + '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.): + """ + Constructor of the Scattered_light_disk object, taking in input the + geometric parameters of the disk, the radial density distribution + and the scattering phase function. + 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 + 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). + + Parameters + ---------- + nx : int + number of pixels along the x axis of the image (default 200) + ny : int + number of pixels along the y axis of the image (default 200) + distance : float + distance to the star in pc (default 70.) + itilt : float + inclination wrt the line of sight in degrees (0 means pole-on, + 90 means edge-on, default 60 degrees) + omega : float + argument of the pericenter in degrees (0 by default) + pxInArcsec : float + pixel field of view in arcsec/px (default the SPHERE pixel + scale 0.01225 arcsec/px) + pa : float + position angle of the disc in degrees (default 0 degrees, e.g. North) + flux_max : float + the max flux of the disk in ADU. By default None, meaning that + the disk flux is not normalized to any value. + density_dico : dict + Parameters describing the dust density distribution function + to be implemented. By default, it uses a two-power law dust + distribution with a vertical gaussian distribution with + linear flaring. This dictionary should at least contain the key + "name". + For a to-power law distribution, you can set it with + 'name:'2PowerLaws' and with the following parameters: + a : float + reference radius in au (default 40) + ksi0 : float + scale height in au at the reference radius (default 1 a.u.) + gamma : float + exponent (2=gaussian,1=exponential profile, default 2) + beta : float + flaring index (0=no flaring, 1=linear flaring, default 1) + ain : float + slope of the power-low distribution in the inner disk. It + must be positive (default 5) + aout : float + slope of the power-low distribution in the outer disk. It + must be negative (default -5) + e : float + eccentricity (default 0) + amin: float + minimim semi-major axis: the dust density is 0 below this + value (default 0) + 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". + xdo : float + disk offset along the x-axis in the disk frame (=semi-major axis), + in a.u. (default 0) + ydo : float + disk offset along the y-axis in the disk frame (=semi-minor axis), + in a.u. (default 0) + """ + 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 + self.distance = distance # distance to the star in pc + self.set_inclination(itilt) + self.set_omega(omega) + self.set_flux_max(flux_max) + self.pxInArcsec = pxInArcsec # pixel field of view in arcsec/px + self.pxInAU = self.pxInArcsec*self.distance # 1 pixel in AU + # disk offset along the x-axis in the disk frame (semi-major axis), AU + self.xdo = xdo + # disk offset along the y-axis in the disk frame (semi-minor axis), AU + self.ydo = ydo + 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))) + 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, + self.y_vector) + self.set_pa(pa) + self.phase_function = Phase_function(spf_dico=spf_dico) + self.scattered_light_map = np.ndarray((ny, nx)) + self.scattered_light_map.fill(0.) + + def set_inclination(self, itilt): + """ + Sets the inclination of the disk. + + Parameters + ---------- + itilt : float + inclination of the disk wrt the line of sight in degrees (0 means + pole-on, 90 means edge-on, default 60 degrees) + """ + self.itilt = float(itilt) # inclination wrt the line of sight in deg + self.cosi = np.cos(np.deg2rad(self.itilt)) + self.sini = np.sin(np.deg2rad(self.itilt)) + + def set_pa(self, pa): + """ + Sets the disk position angle + + Parameters + ---------- + pa : float + position angle in degrees + """ + self.pa = pa # position angle of the disc in degrees + self.cospa = np.cos(np.deg2rad(self.pa)) + self.sinpa = np.sin(np.deg2rad(self.pa)) + # rotation to get the disk major axis properly oriented, x in AU + self.y_map = (self.cospa*self.x_map_0PA + self.sinpa*self.y_map_0PA) + # rotation to get the disk major axis properly oriented, y in AU + self.x_map = (-self.sinpa*self.x_map_0PA + self.cospa*self.y_map_0PA) + + def set_omega(self, omega): + """ + Sets the argument of pericenter + + Parameters + ---------- + omega : float + angle in degrees + """ + self.omega = float(omega) + + def set_flux_max(self, flux_max): + """ + Sets the mas flux of the disk + + Parameters + ---------- + flux_max : float + the max flux of the disk in ADU + """ + self.flux_max = flux_max + + def set_density_distribution(self, density_dico): + """ + Sets or updates the parameters of the density distribution + + Parameters + ---------- + density_dico : dict + Parameters describing the dust density distribution function + to be implemented. By default, it uses a two-power law dust + distribution with a vertical gaussian distribution with + linear flaring. This dictionary should at least contain the key + "name". For a to-power law distribution, you can set it with + name:'2PowerLaws' and with the following parameters: + a : float + reference radius in au (default 60) + ksi0 : float + scale height in au at the reference radius (default 1 a.u.) + gamma : float + exponent (2=gaussian,1=exponential profile, default 2) + beta : float + flaring index (0=no flaring, 1=linear flaring, default 1) + ain : float + slope of the power-low distribution in the inner disk. It + must be positive (default 5) + aout : float + slope of the power-low distribution in the outer disk. It + must be negative (default -5) + e : float + eccentricity (default 0) + """ + self.dust_density.set_density_distribution(density_dico) + + def set_phase_function(self, spf_dico): + """ + Sets the phase function of the dust + + Parameters + ---------- + spf_dico : dict + Parameters describing the scattering phase function to be + implemented. Three phase functions are implemented so far: single + Heyney Greenstein, double Heyney Greenstein and custum phase + functions through interpolation. Read the constructor of each of + those classes to know which parameters must be set in the dictionary + in each case. + """ + self.phase_function = Phase_function(spf_dico=spf_dico) + + def print_info(self): + """ + Prints the information of the disk and image parameters + """ + print('-----------------------------------') + print('Geometrical properties of the image') + print('-----------------------------------') + print('Image size: {0:d} px by {1:d} px'.format(self.nx,self.ny)) + msg1 = 'Pixel size: {0:.4f} arcsec/px or {1:.2f} au/px' + print(msg1.format(self.pxInArcsec, self.pxInAU)) + msg2 = 'Distance of the star {0:.1f} pc' + print(msg2.format(self.distance)) + msg3 = 'From {0:.1f} au to {1:.1f} au in X' + print(msg3.format(self.x_vector[0],self.x_vector[self.nx-1])) + msg4 = 'From {0:.1f} au to {1:.1f} au in Y' + print(msg4.format(self.y_vector[0], self.y_vector[self.nx-1])) + print('Position angle of the disc: {0:.2f} degrees'.format(self.pa)) + print('Inclination {0:.2f} degrees'.format(self.itilt)) + print('Argument of pericenter {0:.2f} degrees'.format(self.omega)) + if self.flux_max is not None: + print('Maximum flux of the disk {0:.2f}'.format(self.flux_max)) + self.dust_density.print_info() + self.phase_function.print_info() + + def check_inclination(self): + """ + Checks whether the inclination set is close to edge-on and risks to + induce artefacts from the limited numerical accuracy. In such a case + the inclination is changed to be less edge-on. + """ + if np.abs(np.mod(self.itilt,180)-90) < np.abs(np.mod(self.dust_density.dust_distribution_calc.itiltthreshold,180)-90): + print('Warning the disk is too close to edge-on') + msg = 'The inclination was changed from {0:.2f} to {1:.2f}' + print(msg.format(self.itilt, + self.dust_density.dust_distribution_calc.itiltthreshold)) + self.set_inclination(self.dust_density.dust_distribution_calc.itiltthreshold) + + def compute_scattered_light(self, halfNbSlices=25): + """ + Computes the scattered lignt image of the disk. + + Parameters + ---------- + halfNbSlices : integer + half number of distances along the line of sight l + """ + self.check_inclination() + # 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)) + # dist to reach +zmax, AU: + lzp_map = self.dust_density.dust_distribution_calc.zmax/self.cosi + \ + lz0_map + # dist to reach -zmax, AU: + lzm_map = -self.dust_density.dust_distribution_calc.zmax/self.cosi + \ + lz0_map + dl_map = np.absolute(lzp_map-lzm_map) # l range, in AU + # squared maximum l value to reach the outer disk radius, in AU^2: + lmax2 = self.dust_density.dust_distribution_calc.rmax**2 - (self.x_map**2+self.y_map**2) + # squared minimum l value to reach the inner disk radius, in AU^2: + lmin2 = (self.x_map**2+self.y_map**2)-self.rmin**2 + validPixel_map = (lmax2 > 0.) * (lmin2 > 0.) + lwidth = 100. # control the distribution of distances along l + nbSlices = 2*halfNbSlices-1 # total number of distances + # along the line of sight + tmp = (np.exp(np.arange(halfNbSlices)*np.log(lwidth+1.) / + (halfNbSlices-1.))-1.)/lwidth # between 0 and 1 + ll = np.concatenate((-tmp[:0:-1], tmp)) + # 1d array pre-calculated values, AU + ycs_vector = self.cosi*self.y_map[validPixel_map] + # 1d array pre-calculated values, AU + zsn_vector = -self.sini*self.y_map[validPixel_map] + xd_vector = self.x_map[validPixel_map] # x_disk, in AU + limage = np.ndarray((nbSlices, self.ny, self.nx)) + limage.fill(0.) + + for il in range(nbSlices): + # distance along the line of sight to reach the plane z + l_vector =lz0_map[validPixel_map] + ll[il]*dl_map[validPixel_map] + # rotation about x axis + yd_vector = ycs_vector + self.sini * l_vector # y_Disk in AU + zd_vector = zsn_vector + self.cosi * l_vector # 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 + # Polar coordinates in the disk frame, and semi-major axis: + # midplane distance to the disk center (r coordinate), in AU + r_vector = np.sqrt((xd_vector-self.xdo)**2+(yd_vector-self.ydo)**2) + # polar angle in radians between 0 and pi + theta_vector = np.arctan2(yd_vector-self.ydo,xd_vector-self.xdo) + costheta_vector = np.cos(theta_vector-np.deg2rad(self.omega)) + # Scattered light: + # volume density + rho_vector = self.dust_density.density_cylindrical(r_vector, + costheta_vector, + zd_vector) + phase_function = self.phase_function.compute_phase_function_from_cosphi(cosphi_vector) + image = np.ndarray((self.ny, self.nx)) + image.fill(0.) + image[validPixel_map] = rho_vector*phase_function/d2star_vector + 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,:,:]) + if self.flux_max is not None: + self.scattered_light_map *= (self.flux_max/np.nanmax(self.scattered_light_map)) + return self.scattered_light_map + + + +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.,'amin':0.,'dens_at_r0':1.}): + """ + Constructor for the Dust_distribution class. + + We assume the dust density is 0 radially after it drops below 0.5% + (the accuracy variable) of the peak density in + 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): + errmsg = 'The parameters describing the dust density distribution' \ + ' must be a Python dictionnary' + raise TypeError(errmsg) + if 'name' not in density_dico.keys(): + 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) + else: + errmsg = 'The only dust distribution implemented so far is the' \ + ' "2PowerLaws"' + raise TypeError(errmsg) + + def set_density_distribution(self,density_dico): + """ + Update the parameters of the density distribution. + """ + self.dust_distribution_calc.set_density_distribution(density_dico) + + def density_cylindrical(self, r, costheta, z): + """ + Return the particule volume density at r, theta, z. + """ + return self.dust_distribution_calc.density_cylindrical(r, costheta, z) + + def density_cartesian(self, x, y, z): + """ + 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) + + def print_info(self, pxInAu=None): + """ + 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) + + +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.,\ + 'amin':0.,'dens_at_r0':1.}): + """ + Constructor for the Dust_distribution class. + + We assume the dust density is 0 radially after it drops below 0.5% + (the accuracy variable) of the peak density in + the midplane, and vertically whenever it drops below 0.5% of the + peak density in the midplane + """ + self.accuracy = accuracy + self.set_density_distribution(density_dico) + + def set_density_distribution(self,density_dico): + """ + """ + if 'ksi0' not in density_dico.keys(): + ksi0 = 1. + else: + ksi0 = density_dico['ksi0'] + if 'beta' not in density_dico.keys(): + beta = 1. + else: + beta = density_dico['beta'] + if 'gamma' not in density_dico.keys(): + gamma = 1. + else: + gamma = density_dico['gamma'] + if 'aout' not in density_dico.keys(): + aout = -5. + else: + aout = density_dico['aout'] + if 'ain' not in density_dico.keys(): + ain = 5. + else: + ain = density_dico['ain'] + if 'e' not in density_dico.keys(): + e = 0. + else: + e = density_dico['e'] + if 'a' not in density_dico.keys(): + a = 60. + else: + a = density_dico['a'] + if 'amin' not in density_dico.keys(): + amin = 0. + else: + amin = density_dico['amin'] + if 'dens_at_r0' not in density_dico.keys(): + dens_at_r0=1. + else: + dens_at_r0=density_dico['dens_at_r0'] + self.set_vertical_density(ksi0=ksi0, gamma=gamma, beta=beta) + self.set_radial_density(ain=ain, aout=aout, a=a, e=e,amin=amin,dens_at_r0=dens_at_r0) + + def set_vertical_density(self, ksi0=1., gamma=2., beta=1.): + """ + Sets the parameters of the vertical density function + + Parameters + ---------- + ksi0 : float + scale height in au at the reference radius (default 1 a.u.) + gamma : float + exponent (2=gaussian,1=exponential profile, default 2) + beta : float + flaring index (0=no flaring, 1=linear flaring, default 1) + """ + if gamma < 0.: + print('Warning the vertical exponent gamma is negative') + print('Gamma was changed from {0:6.2f} to 0.1'.format(gamma)) + gamma = 0.1 + if ksi0 < 0.: + print('Warning the scale height ksi0 is negative') + print('ksi0 was changed from {0:6.2f} to 0.1'.format(ksi0)) + ksi0 = 0.1 + if beta < 0.: + 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) + + def set_radial_density(self, ain=5., aout=-5., a=60., e=0.,amin=0.,dens_at_r0=1.): + """ + Sets the parameters of the radial density function + + Parameters + ---------- + ain : float + slope of the power-low distribution in the inner disk. It + must be positive (default 5) + aout : float + slope of the power-low distribution in the outer disk. It + must be negative (default -5) + a : float + reference radius in au (default 60) + e : float + eccentricity (default 0) + amin: float + minimim semi-major axis: the dust density is 0 below this + value (default 0) + """ + if ain < 0.1: + print('Warning the inner slope is greater than 0.1') + print('ain was changed from {0:6.2f} to 0.1'.format(ain)) + ain = 0.1 + if aout > -0.1: + print('Warning the outer slope is greater than -0.1') + print('aout was changed from {0:6.2f} to -0.1'.format(aout)) + aout = -0.1 + if e < 0: + print('Warning the eccentricity is negative') + print('e was changed from {0:6.2f} to 0'.format(e)) + e = 0. + if e >= 1: + print('Warning the eccentricity is greater or equal to 1') + print('e was changed from {0:6.2f} to 0.99'.format(e)) + e = 0.99 + if a < 0: + raise ValueError('Warning the semi-major axis a is negative') + if amin < 0: + raise ValueError('Warning the minimum radius a is negative') + print('amin was changed from {0:6.2f} to 0.'.format(amin)) + amin = 0. + if dens_at_r0 <0: + raise ValueError('Warning the reference dust density at r0 is negative') + print('It was changed from {0:6.2f} to 1.'.format(dens_at_r0)) + dens_at_r0 = 1. + 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.amin = float(amin) + self.pmin = self.amin*(1-self.e**2) ## we assume the inner hole is also elliptic (convention) + self.dens_at_r0 = float(dens_at_r0) + try: + # 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))) + Gamma_in = self.ain+self.beta + Gamma_out = self.aout+self.beta + self.apeak_surface_density = self.a * np.power(-Gamma_in/Gamma_out, + 1./(2.*(Gamma_in-Gamma_out))) + else: + self.apeak = self.a + self.apeak_surface_density = self.a + except OverflowError: + print('The error occured during the calculation of rmax or apeak') + print('Inner slope: {0:.6e}'.format(self.ain)) + print('Outer slope: {0:.6e}'.format(self.aout)) + print('Accuracy: {0:.6e}'.format(self.accuracy)) + raise OverflowError + except ZeroDivisionError: + print('The error occured during the calculation of rmax or apeak') + print('Inner slope: {0:.6e}'.format(self.ain)) + print('Outer slope: {0:.6e}'.format(self.aout)) + print('Accuracy: {0:.6e}'.format(self.accuracy)) + raise ZeroDivisionError + self.itiltthreshold = np.rad2deg(np.arctan(self.rmax/self.zmax)) + + def print_info(self, pxInAu=None): + """ + Utility function that displays the parameters of the radial distribution + of the dust + + Input: + - pxInAu (optional): the pixel size in au + """ + def rad_density(r): + return np.sqrt(2/(np.power(r/self.a,-2*self.ain) + + np.power(r/self.a,-2*self.aout))) + half_max_density = lambda r:rad_density(r)/rad_density(self.apeak)-1./2. + try: + if self.aout < -3: + a_plus_hwhm = newton(half_max_density,self.apeak*1.04) + else: + a_plus_hwhm = newton(half_max_density,self.apeak*1.1) + except RuntimeError: + a_plus_hwhm = np.nan + try: + if self.ain < 2: + a_minus_hwhm = newton(half_max_density,self.apeak*0.5) + else: + a_minus_hwhm = newton(half_max_density,self.apeak*0.95) + except RuntimeError: + a_minus_hwhm = np.nan + if pxInAu is not None: + 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 in plane z=0: {0:.1f}au or ' \ + '{1:.1f}px (same as ref sma if ain=-aout)' + print(msg2.format(self.apeak,self.apeak/pxInAu)) + msg3 = 'Semi-major axis at half max dust density in plane z=0: {0:.1f}au or ' \ + '{1:.1f}px for the inner edge ' \ + '/ {2:.1f}au or {3:.1f}px for the outer edge, with a FWHM of ' \ + '{4:.1f}au or {5:.1f}px' + print(msg3.format(a_minus_hwhm,a_minus_hwhm/pxInAu,a_plus_hwhm,\ + a_plus_hwhm/pxInAu,a_plus_hwhm-a_minus_hwhm,\ + (a_plus_hwhm-a_minus_hwhm)/pxInAu)) + msg4 = 'Semi-major axis at maximum dust surface density: {0:.1f}au or ' \ + '{1:.1f}px (same as ref sma if ain=-aout)' + print(msg4.format(self.apeak_surface_density,self.apeak_surface_density/pxInAu)) + msg5 = 'Ellipse p parameter: {0:.1f}au or {1:.1f}px' + print(msg5.format(self.p,self.p/pxInAu)) + else: + print('Reference semi-major axis: {0:.1f}au'.format(self.a)) + msg = 'Semi-major axis at maximum dust density in plane z=0: {0:.1f}au (same ' \ + 'as ref sma if ain=-aout)' + print(msg.format(self.apeak)) + msg3 = 'Semi-major axis at half max dust density: {0:.1f}au ' \ + '/ {1:.1f}au for the inner/outer edge, or a FWHM of ' \ + '{2:.1f}au' + print(msg3.format(a_minus_hwhm,a_plus_hwhm,a_plus_hwhm-a_minus_hwhm)) + 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)) + print('Density at the reference semi-major axis: {0:4.3e} (arbitrary unit'.format(self.dens_at_r0)) + if self.amin>0: + print('Minimum radius (sma): {0:.2f}au'.format(self.amin)) + if pxInAu is not None: + 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('Vertical profile index: {0:.2f}'.format(self.gamma)) + 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') + print('------------------------------------') + 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)) + msg = 'Inclination threshold: {0:.2f} degrees' + print(msg.format(self.itiltthreshold)) + return + + 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)) + den = (np.power(radial_ratio, -2*self.ain) + + np.power(radial_ratio,-2*self.aout)) + radial_density_term = np.sqrt(2./den)*self.dens_at_r0 + if self.pmin>0: + radial_density_term[r/(self.pmin/(1-self.e*costheta)) <= 1]=0 + 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 + """ + r = np.sqrt(x**2+y**2) + if r == 0: + costheta=0 + else: + costheta = x/r + return self.density_cylindrical(r,costheta,z) + + +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}): + """ + Constructor of the Phase_function class. It checks whether the spf_dico + contains a correct name and sets the attribute phase_function_calc + + 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 + '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 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 ' \ + 'Python dictionnary' + raise TypeError(msg) + if 'name' not in spf_dico.keys(): + 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 + 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': + self.phase_function_calc = DoubleHenyeyGreenstein_SPF(spf_dico) + elif self.type == 'interpolated': + self.phase_function_calc = Interpolated_SPF(spf_dico) + else: + msg = 'Type of phase function not understood: {0:s}' + raise TypeError(msg.format(self.type)) + + 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. + + Parameters + ---------- + cos_phi : float or array + 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: + 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 + + def print_info(self): + """ + Prints information on the type and parameters of the scattering phase + function. + """ + print('----------------------------') + print('Phase function parameters') + print('----------------------------') + print('Type of phase function: {0:s}'.format(self.type)) + print('Linear polarisation: {0!r}'.format(self.polar)) + self.phase_function_calc.print_info() + + def plot_phase_function(self): + """ + Plots the scattering phase function + """ + phi = np.arange(0, 180, 1) + phase_func = self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi))) + if self.polar: + 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) + plt.plot(phi, phase_func) + plt.xlabel('Scattering phase angle in degrees') + plt.ylabel('Scattering phase function') + plt.grid() + plt.xlim(0, 180) + plt.show() + + +class HenyeyGreenstein_SPF(object): + """ + Implementation of a scattering phase function with a single Henyey + Greenstein function. + """ + + def __init__(self, spf_dico={'g':0.}): + """ + Constructor of a Heyney Greenstein phase function. + + Parameters + ---------- + spf_dico : dictionnary containing the key "g" (float) + g is the Heyney Greenstein coefficient and should be between -1 + (backward scattering) and 1 (forward scattering). + """ + # 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"') + # 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') + self.set_phase_function(spf_dico['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('The value was changed from {0:6.2f} to 0.99'.format(g)) + g = 0.99 + elif g <= -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) + + 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. + + Parameters + ---------- + cos_phi : float or array + 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.) + + def print_info(self): + """ + Prints the value of the HG coefficient + """ + print('Heynyey Greenstein coefficient: {0:.2f}'.format(self.g)) + + +class DoubleHenyeyGreenstein_SPF(object): + """ + 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}): + """ + """ + # 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"') + # 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') + # 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"') + # 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 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') + if len(spf_dico['g']) != 2: + raise TypeError('The keys "weight" and "g" must contain the same' + ' number of elements') + self.g = spf_dico['g'] + self.weight = spf_dico['weight'] + + 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.)) + + 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. + + 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. + """ + 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. + + Parameters + ---------- + cos_phi : float or array + 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) + + +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])}): + """ + Constructor of the Interpolated_SPF class. It checks whether the spf_dico + contains the keys 'phi' and 'spf' + + Parameters + ---------- + spf_dico : dict + dictionnary containing at least the keys "phi" (list of scattering angles) + and "spf" (list of corresponding scattering phase function values) + Optionnaly it can specify the kind of interpolation requested (as + specified by the scipy.interpolate.interp1d function), by default + it uses a quadratic interpolation. + """ + 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)) + elif not isinstance(spf_dico[key],(list,tuple,np.ndarray)): + raise TypeError('The key "{0:s}" of a "interpolated" phase' + ' function dictionnary must be a list, np 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') + self.interpolate_phase_function(spf_dico) + + def print_info(self): + """ + Prints the information of the spf + """ + 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) + + 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. By default it uses a quadractic interpolation. + + Parameters + ---------- + spf_dico : dict + dictionnary containing at least the keys "phi" (list of scattering angles) + and "spf" (list of corresponding scattering phase function values) + Optionnaly it can specify the kind of interpolation requested (as + specified by the scipy.interpolate.interp1d function), by default + it uses a quadratic interpolation. + """ + if 'kind' in spf_dico.keys(): + if not isinstance(spf_dico['kind'],int) and spf_dico['kind'] not in \ + ['linear', 'nearest', 'zero', 'slinear', \ + 'quadratic', 'cubic', 'previous', 'next']: + 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'])) + 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) + # 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): + """ + 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. + """ + return self.interpolation_function(np.rad2deg(np.arccos(cos_phi))) + # return self.interpolation_function(cos_phi) + + +if __name__ == '__main__': + """ + It is an example of use of the ScatteredLightDisk class. + """ + + # Example 1: creation of a disk at 20pc, with semi-major axis 20 a.u. + # inner and outer slopes 2 and -5, inclined by 70degrees wrt + # the line of sight and with a PA of 45degrees, an isotropic + # phase function + Scattered_light_disk1 = ScatteredLightDisk(nx=201, ny=201, distance=20, + itilt=70.,omega=0., + pxInArcsec=0.01225, pa=45., + density_dico={'name':'2PowerLaws', + 'ain':5,'aout':-5, + 'a':40,'e':0, + 'ksi0':1., + 'gamma':2., + 'beta':1.}, + spf_dico={'name':'HG','g':0., + 'polar':False}) + disk1 = Scattered_light_disk1.compute_scattered_light() + + # If you prefer set the parameters in differen steps, you can also do that + # with: (this is identical) + Scattered_light_disk1 = ScatteredLightDisk(nx=201, ny=201, distance=20) + Scattered_light_disk1.set_pa(45) + Scattered_light_disk1.set_inclination(70) + Scattered_light_disk1.set_density_distribution({'name':'2PowerLaws', + 'ain':2,'aout':-5,'a':20, + 'e':0,'ksi0':1., + 'gamma':2.,'beta':1.}) + Scattered_light_disk1.set_phase_function({'name':'HG','g':0.,'polar':False}) + disk1 = Scattered_light_disk1.compute_scattered_light() + + # If you want to know what are the parameters you set for your disk: + Scattered_light_disk1.print_info() + + # Example 2: Creation of a disk similar to example 1 but with a double + # Heyney Greenstein phase function + + Scattered_light_disk2 = ScatteredLightDisk(nx=201,ny=201,distance=20) + Scattered_light_disk2.set_pa(45) + Scattered_light_disk2.set_inclination(70) + Scattered_light_disk2.set_density_distribution({'name':'2PowerLaws', + 'ain':2,'aout':-5,'a':20, + 'e':0,'ksi0':1., + 'gamma':2.,'beta':1.}) + Scattered_light_disk2.set_phase_function({'name':'DoubleHG', + 'g':[0.6,-0.6],'weight':0.7, + 'polar':False}) + Scattered_light_disk2.print_info() + disk2 = Scattered_light_disk2.compute_scattered_light() + + # Let's turn the polarisation on now: + Scattered_light_disk2.set_phase_function({'name':'DoubleHG', + 'g':[0.6,-0.6],'weight':0.7, + 'polar':True}) + disk2_polar = Scattered_light_disk2.compute_scattered_light() + diff --git a/vip_hci/negfc/utils_mcmc.py b/vip_hci/fm/utils_mcmc.py similarity index 100% rename from vip_hci/negfc/utils_mcmc.py rename to vip_hci/fm/utils_mcmc.py diff --git a/vip_hci/negfc/utils_negfc.py b/vip_hci/fm/utils_negfc.py similarity index 98% rename from vip_hci/negfc/utils_negfc.py rename to vip_hci/fm/utils_negfc.py index 8077861d..6f7d145c 100644 --- a/vip_hci/negfc/utils_negfc.py +++ b/vip_hci/fm/utils_negfc.py @@ -9,7 +9,7 @@ __all__ = ['cube_planet_free'] import numpy as np -from ..metrics import cube_inject_companions +from ..fm import cube_inject_companions def cube_planet_free(planet_parameter, cube, angs, psfn, plsc, imlib='vip-fft', diff --git a/vip_hci/hci_dataset.py b/vip_hci/hci_dataset.py index 18fb7c49..05d91991 100644 --- a/vip_hci/hci_dataset.py +++ b/vip_hci/hci_dataset.py @@ -12,6 +12,8 @@ import copy import hciplot as hp from .fits import open_fits +from .fm import (cube_inject_companions, generate_cube_copies_with_injections, + normalize_psf) from .preproc import (frame_crop, frame_px_resampling, frame_rotate, frame_shift, frame_center_satspots, frame_center_radon) from .preproc import (cube_collapse, cube_crop_frames, cube_derotate, @@ -24,9 +26,7 @@ cube_filter_highpass, cube_filter_lowpass, mask_circle) from .stats import (frame_basic_stats, frame_histo_stats, frame_average_radprofile, cube_basic_stats, cube_distance) -from .metrics import (frame_report, cube_inject_companions, - generate_cube_copies_with_injections, snr, - snrmap, detection, normalize_psf) +from .metrics import (frame_report, snr, snrmap, detection) from .config.utils_conf import check_array, Saveable, print_precision from .config.mem import check_enough_memory diff --git a/vip_hci/metrics/__init__.py b/vip_hci/metrics/__init__.py index a9a4d54d..3e5ca512 100644 --- a/vip_hci/metrics/__init__.py +++ b/vip_hci/metrics/__init__.py @@ -12,9 +12,6 @@ """ from .contrcurve import * from .detection import * -from .fakecomp import * from .roc import * from .snr_source import * -from .fakedisk import * -from .scattered_light_disk import * from .stim import * diff --git a/vip_hci/metrics/contrcurve.py b/vip_hci/metrics/contrcurve.py index 1c492ac4..6c652a4c 100644 --- a/vip_hci/metrics/contrcurve.py +++ b/vip_hci/metrics/contrcurve.py @@ -19,8 +19,8 @@ from scipy.signal import savgol_filter from skimage.draw import disk from matplotlib import pyplot as plt -from .fakecomp import (cube_inject_companions, frame_inject_companion, - normalize_psf) +from ..fm import (cube_inject_companions, frame_inject_companion, + normalize_psf) from ..config import time_ini, timing from ..config.utils_conf import sep from ..var import frame_center, dist diff --git a/vip_hci/metrics/dust_distribution.py b/vip_hci/metrics/dust_distribution.py deleted file mode 100644 index f25b71d8..00000000 --- a/vip_hci/metrics/dust_distribution.py +++ /dev/null @@ -1,365 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed May 6 17:07:00 2015 -""" -import numpy as np -from scipy.optimize import newton - -_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.,'amin':0.,'dens_at_r0':1.}): - """ - Constructor for the Dust_distribution class. - - We assume the dust density is 0 radially after it drops below 0.5% - (the accuracy variable) of the peak density in - 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): - errmsg = 'The parameters describing the dust density distribution' \ - ' must be a Python dictionnary' - raise TypeError(errmsg) - if 'name' not in density_dico.keys(): - 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) - else: - errmsg = 'The only dust distribution implemented so far is the' \ - ' "2PowerLaws"' - raise TypeError(errmsg) - - def set_density_distribution(self,density_dico): - """ - Update the parameters of the density distribution. - """ - self.dust_distribution_calc.set_density_distribution(density_dico) - - def density_cylindrical(self, r, costheta, z): - """ - Return the particule volume density at r, theta, z. - """ - return self.dust_distribution_calc.density_cylindrical(r, costheta, z) - - def density_cartesian(self, x, y, z): - """ - 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) - - def print_info(self, pxInAu=None): - """ - 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) - - -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.,\ - 'amin':0.,'dens_at_r0':1.}): - """ - Constructor for the Dust_distribution class. - - We assume the dust density is 0 radially after it drops below 0.5% - (the accuracy variable) of the peak density in - the midplane, and vertically whenever it drops below 0.5% of the - peak density in the midplane - """ - self.accuracy = accuracy - self.set_density_distribution(density_dico) - - def set_density_distribution(self,density_dico): - """ - """ - if 'ksi0' not in density_dico.keys(): - ksi0 = 1. - else: - ksi0 = density_dico['ksi0'] - if 'beta' not in density_dico.keys(): - beta = 1. - else: - beta = density_dico['beta'] - if 'gamma' not in density_dico.keys(): - gamma = 1. - else: - gamma = density_dico['gamma'] - if 'aout' not in density_dico.keys(): - aout = -5. - else: - aout = density_dico['aout'] - if 'ain' not in density_dico.keys(): - ain = 5. - else: - ain = density_dico['ain'] - if 'e' not in density_dico.keys(): - e = 0. - else: - e = density_dico['e'] - if 'a' not in density_dico.keys(): - a = 60. - else: - a = density_dico['a'] - if 'amin' not in density_dico.keys(): - amin = 0. - else: - amin = density_dico['amin'] - if 'dens_at_r0' not in density_dico.keys(): - dens_at_r0=1. - else: - dens_at_r0=density_dico['dens_at_r0'] - self.set_vertical_density(ksi0=ksi0, gamma=gamma, beta=beta) - self.set_radial_density(ain=ain, aout=aout, a=a, e=e,amin=amin,dens_at_r0=dens_at_r0) - - def set_vertical_density(self, ksi0=1., gamma=2., beta=1.): - """ - Sets the parameters of the vertical density function - - Parameters - ---------- - ksi0 : float - scale height in au at the reference radius (default 1 a.u.) - gamma : float - exponent (2=gaussian,1=exponential profile, default 2) - beta : float - flaring index (0=no flaring, 1=linear flaring, default 1) - """ - if gamma < 0.: - print('Warning the vertical exponent gamma is negative') - print('Gamma was changed from {0:6.2f} to 0.1'.format(gamma)) - gamma = 0.1 - if ksi0 < 0.: - print('Warning the scale height ksi0 is negative') - print('ksi0 was changed from {0:6.2f} to 0.1'.format(ksi0)) - ksi0 = 0.1 - if beta < 0.: - 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) - - def set_radial_density(self, ain=5., aout=-5., a=60., e=0.,amin=0.,dens_at_r0=1.): - """ - Sets the parameters of the radial density function - - Parameters - ---------- - ain : float - slope of the power-low distribution in the inner disk. It - must be positive (default 5) - aout : float - slope of the power-low distribution in the outer disk. It - must be negative (default -5) - a : float - reference radius in au (default 60) - e : float - eccentricity (default 0) - amin: float - minimim semi-major axis: the dust density is 0 below this - value (default 0) - """ - if ain < 0.1: - print('Warning the inner slope is greater than 0.1') - print('ain was changed from {0:6.2f} to 0.1'.format(ain)) - ain = 0.1 - if aout > -0.1: - print('Warning the outer slope is greater than -0.1') - print('aout was changed from {0:6.2f} to -0.1'.format(aout)) - aout = -0.1 - if e < 0: - print('Warning the eccentricity is negative') - print('e was changed from {0:6.2f} to 0'.format(e)) - e = 0. - if e >= 1: - print('Warning the eccentricity is greater or equal to 1') - print('e was changed from {0:6.2f} to 0.99'.format(e)) - e = 0.99 - if a < 0: - raise ValueError('Warning the semi-major axis a is negative') - if amin < 0: - raise ValueError('Warning the minimum radius a is negative') - print('amin was changed from {0:6.2f} to 0.'.format(amin)) - amin = 0. - if dens_at_r0 <0: - raise ValueError('Warning the reference dust density at r0 is negative') - print('It was changed from {0:6.2f} to 1.'.format(dens_at_r0)) - dens_at_r0 = 1. - 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.amin = float(amin) - self.pmin = self.amin*(1-self.e**2) ## we assume the inner hole is also elliptic (convention) - self.dens_at_r0 = float(dens_at_r0) - try: - # 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))) - Gamma_in = self.ain+self.beta - Gamma_out = self.aout+self.beta - self.apeak_surface_density = self.a * np.power(-Gamma_in/Gamma_out, - 1./(2.*(Gamma_in-Gamma_out))) - else: - self.apeak = self.a - self.apeak_surface_density = self.a - except OverflowError: - print('The error occured during the calculation of rmax or apeak') - print('Inner slope: {0:.6e}'.format(self.ain)) - print('Outer slope: {0:.6e}'.format(self.aout)) - print('Accuracy: {0:.6e}'.format(self.accuracy)) - raise OverflowError - except ZeroDivisionError: - print('The error occured during the calculation of rmax or apeak') - print('Inner slope: {0:.6e}'.format(self.ain)) - print('Outer slope: {0:.6e}'.format(self.aout)) - print('Accuracy: {0:.6e}'.format(self.accuracy)) - raise ZeroDivisionError - self.itiltthreshold = np.rad2deg(np.arctan(self.rmax/self.zmax)) - - def print_info(self, pxInAu=None): - """ - Utility function that displays the parameters of the radial distribution - of the dust - - Input: - - pxInAu (optional): the pixel size in au - """ - def rad_density(r): - return np.sqrt(2/(np.power(r/self.a,-2*self.ain) + - np.power(r/self.a,-2*self.aout))) - half_max_density = lambda r:rad_density(r)/rad_density(self.apeak)-1./2. - try: - if self.aout < -3: - a_plus_hwhm = newton(half_max_density,self.apeak*1.04) - else: - a_plus_hwhm = newton(half_max_density,self.apeak*1.1) - except RuntimeError: - a_plus_hwhm = np.nan - try: - if self.ain < 2: - a_minus_hwhm = newton(half_max_density,self.apeak*0.5) - else: - a_minus_hwhm = newton(half_max_density,self.apeak*0.95) - except RuntimeError: - a_minus_hwhm = np.nan - if pxInAu is not None: - 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 in plane z=0: {0:.1f}au or ' \ - '{1:.1f}px (same as ref sma if ain=-aout)' - print(msg2.format(self.apeak,self.apeak/pxInAu)) - msg3 = 'Semi-major axis at half max dust density in plane z=0: {0:.1f}au or ' \ - '{1:.1f}px for the inner edge ' \ - '/ {2:.1f}au or {3:.1f}px for the outer edge, with a FWHM of ' \ - '{4:.1f}au or {5:.1f}px' - print(msg3.format(a_minus_hwhm,a_minus_hwhm/pxInAu,a_plus_hwhm,\ - a_plus_hwhm/pxInAu,a_plus_hwhm-a_minus_hwhm,\ - (a_plus_hwhm-a_minus_hwhm)/pxInAu)) - msg4 = 'Semi-major axis at maximum dust surface density: {0:.1f}au or ' \ - '{1:.1f}px (same as ref sma if ain=-aout)' - print(msg4.format(self.apeak_surface_density,self.apeak_surface_density/pxInAu)) - msg5 = 'Ellipse p parameter: {0:.1f}au or {1:.1f}px' - print(msg5.format(self.p,self.p/pxInAu)) - else: - print('Reference semi-major axis: {0:.1f}au'.format(self.a)) - msg = 'Semi-major axis at maximum dust density in plane z=0: {0:.1f}au (same ' \ - 'as ref sma if ain=-aout)' - print(msg.format(self.apeak)) - msg3 = 'Semi-major axis at half max dust density: {0:.1f}au ' \ - '/ {1:.1f}au for the inner/outer edge, or a FWHM of ' \ - '{2:.1f}au' - print(msg3.format(a_minus_hwhm,a_plus_hwhm,a_plus_hwhm-a_minus_hwhm)) - 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)) - print('Density at the reference semi-major axis: {0:4.3e} (arbitrary unit'.format(self.dens_at_r0)) - if self.amin>0: - print('Minimum radius (sma): {0:.2f}au'.format(self.amin)) - if pxInAu is not None: - 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('Vertical profile index: {0:.2f}'.format(self.gamma)) - 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') - print('------------------------------------') - 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)) - msg = 'Inclination threshold: {0:.2f} degrees' - print(msg.format(self.itiltthreshold)) - return - - 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)) - den = (np.power(radial_ratio, -2*self.ain) + - np.power(radial_ratio,-2*self.aout)) - radial_density_term = np.sqrt(2./den)*self.dens_at_r0 - if self.pmin>0: - radial_density_term[r/(self.pmin/(1-self.e*costheta)) <= 1]=0 - 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 - """ - r = np.sqrt(x**2+y**2) - if r == 0: - costheta=0 - else: - costheta = x/r - return self.density_cylindrical(r,costheta,z) - - -if __name__ == '__main__': - """ - Main of the class for debugging - """ - test = DustEllipticalDistribution2PowerLaws() - test.set_radial_density(ain=5., aout=-5., a=60., e=0.) - test.print_info() - 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 diff --git a/vip_hci/metrics/phase_function.py b/vip_hci/metrics/phase_function.py deleted file mode 100644 index 2f0c1f98..00000000 --- a/vip_hci/metrics/phase_function.py +++ /dev/null @@ -1,378 +0,0 @@ -#! /usr/bin/env python -""" -Phase_function class definition -""" - -__author__ = 'Julien Milli' -__all__ = [] - -import numpy as np -import matplotlib.pyplot as plt -from scipy.interpolate import interp1d - - -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}): - """ - Constructor of the Phase_function class. It checks whether the spf_dico - contains a correct name and sets the attribute phase_function_calc - - 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 - '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 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 ' \ - 'Python dictionnary' - raise TypeError(msg) - if 'name' not in spf_dico.keys(): - 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 - 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': - self.phase_function_calc = DoubleHenyeyGreenstein_SPF(spf_dico) - elif self.type == 'interpolated': - self.phase_function_calc = Interpolated_SPF(spf_dico) - else: - msg = 'Type of phase function not understood: {0:s}' - raise TypeError(msg.format(self.type)) - - 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. - - Parameters - ---------- - cos_phi : float or array - 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: - 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 - - def print_info(self): - """ - Prints information on the type and parameters of the scattering phase - function. - """ - print('----------------------------') - print('Phase function parameters') - print('----------------------------') - print('Type of phase function: {0:s}'.format(self.type)) - print('Linear polarisation: {0!r}'.format(self.polar)) - self.phase_function_calc.print_info() - - def plot_phase_function(self): - """ - Plots the scattering phase function - """ - phi = np.arange(0, 180, 1) - phase_func = self.compute_phase_function_from_cosphi(np.cos(np.deg2rad(phi))) - if self.polar: - 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) - plt.plot(phi, phase_func) - plt.xlabel('Scattering phase angle in degrees') - plt.ylabel('Scattering phase function') - plt.grid() - plt.xlim(0, 180) - plt.show() - - -class HenyeyGreenstein_SPF(object): - """ - Implementation of a scattering phase function with a single Henyey - Greenstein function. - """ - - def __init__(self, spf_dico={'g':0.}): - """ - Constructor of a Heyney Greenstein phase function. - - Parameters - ---------- - spf_dico : dictionnary containing the key "g" (float) - g is the Heyney Greenstein coefficient and should be between -1 - (backward scattering) and 1 (forward scattering). - """ - # 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"') - # 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') - self.set_phase_function(spf_dico['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('The value was changed from {0:6.2f} to 0.99'.format(g)) - g = 0.99 - elif g <= -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) - - 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. - - Parameters - ---------- - cos_phi : float or array - 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.) - - def print_info(self): - """ - Prints the value of the HG coefficient - """ - print('Heynyey Greenstein coefficient: {0:.2f}'.format(self.g)) - - -class DoubleHenyeyGreenstein_SPF(object): - """ - 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}): - """ - """ - # 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"') - # 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') - # 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"') - # 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 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') - if len(spf_dico['g']) != 2: - raise TypeError('The keys "weight" and "g" must contain the same' - ' number of elements') - self.g = spf_dico['g'] - self.weight = spf_dico['weight'] - - 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.)) - - 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. - - 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. - """ - 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. - - Parameters - ---------- - cos_phi : float or array - 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) - - -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])}): - """ - Constructor of the Interpolated_SPF class. It checks whether the spf_dico - contains the keys 'phi' and 'spf' - - Parameters - ---------- - spf_dico : dict - dictionnary containing at least the keys "phi" (list of scattering angles) - and "spf" (list of corresponding scattering phase function values) - Optionnaly it can specify the kind of interpolation requested (as - specified by the scipy.interpolate.interp1d function), by default - it uses a quadratic interpolation. - """ - 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)) - elif not isinstance(spf_dico[key],(list,tuple,np.ndarray)): - raise TypeError('The key "{0:s}" of a "interpolated" phase' - ' function dictionnary must be a list, np 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') - self.interpolate_phase_function(spf_dico) - - def print_info(self): - """ - Prints the information of the spf - """ - 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) - - 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. By default it uses a quadractic interpolation. - - Parameters - ---------- - spf_dico : dict - dictionnary containing at least the keys "phi" (list of scattering angles) - and "spf" (list of corresponding scattering phase function values) - Optionnaly it can specify the kind of interpolation requested (as - specified by the scipy.interpolate.interp1d function), by default - it uses a quadratic interpolation. - """ - if 'kind' in spf_dico.keys(): - if not isinstance(spf_dico['kind'],int) and spf_dico['kind'] not in \ - ['linear', 'nearest', 'zero', 'slinear', \ - 'quadratic', 'cubic', 'previous', 'next']: - 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'])) - 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) - # 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): - """ - 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. - """ - return self.interpolation_function(np.rad2deg(np.arccos(cos_phi))) - # return self.interpolation_function(cos_phi) diff --git a/vip_hci/metrics/roc.py b/vip_hci/metrics/roc.py index 2e6638fc..d98269e4 100644 --- a/vip_hci/metrics/roc.py +++ b/vip_hci/metrics/roc.py @@ -11,12 +11,10 @@ from scipy import stats from photutils import detect_sources from munch import Munch -from ..psfsub.svd import SVDecomposer -from ..var import frame_center, get_annulus_segments from ..config import time_ini, timing, Progressbar -from ..var import get_circle -from .fakecomp import cube_inject_companions - +from ..fm import cube_inject_companions +from ..psfsub.svd import SVDecomposer +from ..var import frame_center, get_annulus_segments, get_circle # TODO: remove the munch dependency class EvalRoc(object): diff --git a/vip_hci/metrics/scattered_light_disk.py b/vip_hci/metrics/scattered_light_disk.py deleted file mode 100644 index 498abdb8..00000000 --- a/vip_hci/metrics/scattered_light_disk.py +++ /dev/null @@ -1,398 +0,0 @@ -#! /usr/bin/env python -""" -Scattered_light_disk class definition -""" - -__author__ = 'Julien Milli' -__all__ = [] - -import numpy as np -from ..var import frame_center -from .dust_distribution import Dust_distribution -from .phase_function import Phase_function - - -class ScatteredLightDisk(object): - """ - Class used to generate a synthetic disc, inspired from a light version of - the GRATER tool (GRenoble RAdiative TransfER, Augereau et al. 1999) - written originally in IDL and converted to Python by J. Milli. - """ - - def __init__(self, nx=200, ny=200, distance=50., itilt=60., omega=0., - pxInArcsec=0.01225, pa=0., flux_max=None, - density_dico={'name': '2PowerLaws', 'ain':5, 'aout':-5, - '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.): - """ - Constructor of the Scattered_light_disk object, taking in input the - geometric parameters of the disk, the radial density distribution - and the scattering phase function. - 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 - 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). - - Parameters - ---------- - nx : int - number of pixels along the x axis of the image (default 200) - ny : int - number of pixels along the y axis of the image (default 200) - distance : float - distance to the star in pc (default 70.) - itilt : float - inclination wrt the line of sight in degrees (0 means pole-on, - 90 means edge-on, default 60 degrees) - omega : float - argument of the pericenter in degrees (0 by default) - pxInArcsec : float - pixel field of view in arcsec/px (default the SPHERE pixel - scale 0.01225 arcsec/px) - pa : float - position angle of the disc in degrees (default 0 degrees, e.g. North) - flux_max : float - the max flux of the disk in ADU. By default None, meaning that - the disk flux is not normalized to any value. - density_dico : dict - Parameters describing the dust density distribution function - to be implemented. By default, it uses a two-power law dust - distribution with a vertical gaussian distribution with - linear flaring. This dictionary should at least contain the key - "name". - For a to-power law distribution, you can set it with - 'name:'2PowerLaws' and with the following parameters: - a : float - reference radius in au (default 40) - ksi0 : float - scale height in au at the reference radius (default 1 a.u.) - gamma : float - exponent (2=gaussian,1=exponential profile, default 2) - beta : float - flaring index (0=no flaring, 1=linear flaring, default 1) - ain : float - slope of the power-low distribution in the inner disk. It - must be positive (default 5) - aout : float - slope of the power-low distribution in the outer disk. It - must be negative (default -5) - e : float - eccentricity (default 0) - amin: float - minimim semi-major axis: the dust density is 0 below this - value (default 0) - 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". - xdo : float - disk offset along the x-axis in the disk frame (=semi-major axis), - in a.u. (default 0) - ydo : float - disk offset along the y-axis in the disk frame (=semi-minor axis), - in a.u. (default 0) - """ - 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 - self.distance = distance # distance to the star in pc - self.set_inclination(itilt) - self.set_omega(omega) - self.set_flux_max(flux_max) - self.pxInArcsec = pxInArcsec # pixel field of view in arcsec/px - self.pxInAU = self.pxInArcsec*self.distance # 1 pixel in AU - # disk offset along the x-axis in the disk frame (semi-major axis), AU - self.xdo = xdo - # disk offset along the y-axis in the disk frame (semi-minor axis), AU - self.ydo = ydo - 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))) - 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, - self.y_vector) - self.set_pa(pa) - self.phase_function = Phase_function(spf_dico=spf_dico) - self.scattered_light_map = np.ndarray((ny, nx)) - self.scattered_light_map.fill(0.) - - def set_inclination(self, itilt): - """ - Sets the inclination of the disk. - - Parameters - ---------- - itilt : float - inclination of the disk wrt the line of sight in degrees (0 means - pole-on, 90 means edge-on, default 60 degrees) - """ - self.itilt = float(itilt) # inclination wrt the line of sight in deg - self.cosi = np.cos(np.deg2rad(self.itilt)) - self.sini = np.sin(np.deg2rad(self.itilt)) - - def set_pa(self, pa): - """ - Sets the disk position angle - - Parameters - ---------- - pa : float - position angle in degrees - """ - self.pa = pa # position angle of the disc in degrees - self.cospa = np.cos(np.deg2rad(self.pa)) - self.sinpa = np.sin(np.deg2rad(self.pa)) - # rotation to get the disk major axis properly oriented, x in AU - self.y_map = (self.cospa*self.x_map_0PA + self.sinpa*self.y_map_0PA) - # rotation to get the disk major axis properly oriented, y in AU - self.x_map = (-self.sinpa*self.x_map_0PA + self.cospa*self.y_map_0PA) - - def set_omega(self, omega): - """ - Sets the argument of pericenter - - Parameters - ---------- - omega : float - angle in degrees - """ - self.omega = float(omega) - - def set_flux_max(self, flux_max): - """ - Sets the mas flux of the disk - - Parameters - ---------- - flux_max : float - the max flux of the disk in ADU - """ - self.flux_max = flux_max - - def set_density_distribution(self, density_dico): - """ - Sets or updates the parameters of the density distribution - - Parameters - ---------- - density_dico : dict - Parameters describing the dust density distribution function - to be implemented. By default, it uses a two-power law dust - distribution with a vertical gaussian distribution with - linear flaring. This dictionary should at least contain the key - "name". For a to-power law distribution, you can set it with - name:'2PowerLaws' and with the following parameters: - a : float - reference radius in au (default 60) - ksi0 : float - scale height in au at the reference radius (default 1 a.u.) - gamma : float - exponent (2=gaussian,1=exponential profile, default 2) - beta : float - flaring index (0=no flaring, 1=linear flaring, default 1) - ain : float - slope of the power-low distribution in the inner disk. It - must be positive (default 5) - aout : float - slope of the power-low distribution in the outer disk. It - must be negative (default -5) - e : float - eccentricity (default 0) - """ - self.dust_density.set_density_distribution(density_dico) - - def set_phase_function(self, spf_dico): - """ - Sets the phase function of the dust - - Parameters - ---------- - spf_dico : dict - Parameters describing the scattering phase function to be - implemented. Three phase functions are implemented so far: single - Heyney Greenstein, double Heyney Greenstein and custum phase - functions through interpolation. Read the constructor of each of - those classes to know which parameters must be set in the dictionary - in each case. - """ - self.phase_function = Phase_function(spf_dico=spf_dico) - - def print_info(self): - """ - Prints the information of the disk and image parameters - """ - print('-----------------------------------') - print('Geometrical properties of the image') - print('-----------------------------------') - print('Image size: {0:d} px by {1:d} px'.format(self.nx,self.ny)) - msg1 = 'Pixel size: {0:.4f} arcsec/px or {1:.2f} au/px' - print(msg1.format(self.pxInArcsec, self.pxInAU)) - msg2 = 'Distance of the star {0:.1f} pc' - print(msg2.format(self.distance)) - msg3 = 'From {0:.1f} au to {1:.1f} au in X' - print(msg3.format(self.x_vector[0],self.x_vector[self.nx-1])) - msg4 = 'From {0:.1f} au to {1:.1f} au in Y' - print(msg4.format(self.y_vector[0], self.y_vector[self.nx-1])) - print('Position angle of the disc: {0:.2f} degrees'.format(self.pa)) - print('Inclination {0:.2f} degrees'.format(self.itilt)) - print('Argument of pericenter {0:.2f} degrees'.format(self.omega)) - if self.flux_max is not None: - print('Maximum flux of the disk {0:.2f}'.format(self.flux_max)) - self.dust_density.print_info() - self.phase_function.print_info() - - def check_inclination(self): - """ - Checks whether the inclination set is close to edge-on and risks to - induce artefacts from the limited numerical accuracy. In such a case - the inclination is changed to be less edge-on. - """ - if np.abs(np.mod(self.itilt,180)-90) < np.abs(np.mod(self.dust_density.dust_distribution_calc.itiltthreshold,180)-90): - print('Warning the disk is too close to edge-on') - msg = 'The inclination was changed from {0:.2f} to {1:.2f}' - print(msg.format(self.itilt, - self.dust_density.dust_distribution_calc.itiltthreshold)) - self.set_inclination(self.dust_density.dust_distribution_calc.itiltthreshold) - - def compute_scattered_light(self, halfNbSlices=25): - """ - Computes the scattered lignt image of the disk. - - Parameters - ---------- - halfNbSlices : integer - half number of distances along the line of sight l - """ - self.check_inclination() - # 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)) - # dist to reach +zmax, AU: - lzp_map = self.dust_density.dust_distribution_calc.zmax/self.cosi + \ - lz0_map - # dist to reach -zmax, AU: - lzm_map = -self.dust_density.dust_distribution_calc.zmax/self.cosi + \ - lz0_map - dl_map = np.absolute(lzp_map-lzm_map) # l range, in AU - # squared maximum l value to reach the outer disk radius, in AU^2: - lmax2 = self.dust_density.dust_distribution_calc.rmax**2 - (self.x_map**2+self.y_map**2) - # squared minimum l value to reach the inner disk radius, in AU^2: - lmin2 = (self.x_map**2+self.y_map**2)-self.rmin**2 - validPixel_map = (lmax2 > 0.) * (lmin2 > 0.) - lwidth = 100. # control the distribution of distances along l - nbSlices = 2*halfNbSlices-1 # total number of distances - # along the line of sight - tmp = (np.exp(np.arange(halfNbSlices)*np.log(lwidth+1.) / - (halfNbSlices-1.))-1.)/lwidth # between 0 and 1 - ll = np.concatenate((-tmp[:0:-1], tmp)) - # 1d array pre-calculated values, AU - ycs_vector = self.cosi*self.y_map[validPixel_map] - # 1d array pre-calculated values, AU - zsn_vector = -self.sini*self.y_map[validPixel_map] - xd_vector = self.x_map[validPixel_map] # x_disk, in AU - limage = np.ndarray((nbSlices, self.ny, self.nx)) - limage.fill(0.) - - for il in range(nbSlices): - # distance along the line of sight to reach the plane z - l_vector =lz0_map[validPixel_map] + ll[il]*dl_map[validPixel_map] - # rotation about x axis - yd_vector = ycs_vector + self.sini * l_vector # y_Disk in AU - zd_vector = zsn_vector + self.cosi * l_vector # 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 - # Polar coordinates in the disk frame, and semi-major axis: - # midplane distance to the disk center (r coordinate), in AU - r_vector = np.sqrt((xd_vector-self.xdo)**2+(yd_vector-self.ydo)**2) - # polar angle in radians between 0 and pi - theta_vector = np.arctan2(yd_vector-self.ydo,xd_vector-self.xdo) - costheta_vector = np.cos(theta_vector-np.deg2rad(self.omega)) - # Scattered light: - # volume density - rho_vector = self.dust_density.density_cylindrical(r_vector, - costheta_vector, - zd_vector) - phase_function = self.phase_function.compute_phase_function_from_cosphi(cosphi_vector) - image = np.ndarray((self.ny, self.nx)) - image.fill(0.) - image[validPixel_map] = rho_vector*phase_function/d2star_vector - 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,:,:]) - if self.flux_max is not None: - self.scattered_light_map *= (self.flux_max/np.nanmax(self.scattered_light_map)) - return self.scattered_light_map - - -if __name__ == '__main__': - """ - It is an example of use of the ScatteredLightDisk class. - """ - - # Example 1: creation of a disk at 20pc, with semi-major axis 20 a.u. - # inner and outer slopes 2 and -5, inclined by 70degrees wrt - # the line of sight and with a PA of 45degrees, an isotropic - # phase function - Scattered_light_disk1 = ScatteredLightDisk(nx=201, ny=201, distance=20, - itilt=70.,omega=0., - pxInArcsec=0.01225, pa=45., - density_dico={'name':'2PowerLaws', - 'ain':5,'aout':-5, - 'a':40,'e':0, - 'ksi0':1., - 'gamma':2., - 'beta':1.}, - spf_dico={'name':'HG','g':0., - 'polar':False}) - disk1 = Scattered_light_disk1.compute_scattered_light() - - # If you prefer set the parameters in differen steps, you can also do that - # with: (this is identical) - Scattered_light_disk1 = ScatteredLightDisk(nx=201, ny=201, distance=20) - Scattered_light_disk1.set_pa(45) - Scattered_light_disk1.set_inclination(70) - Scattered_light_disk1.set_density_distribution({'name':'2PowerLaws', - 'ain':2,'aout':-5,'a':20, - 'e':0,'ksi0':1., - 'gamma':2.,'beta':1.}) - Scattered_light_disk1.set_phase_function({'name':'HG','g':0.,'polar':False}) - disk1 = Scattered_light_disk1.compute_scattered_light() - - # If you want to know what are the parameters you set for your disk: - Scattered_light_disk1.print_info() - - # Example 2: Creation of a disk similar to example 1 but with a double - # Heyney Greenstein phase function - - Scattered_light_disk2 = ScatteredLightDisk(nx=201,ny=201,distance=20) - Scattered_light_disk2.set_pa(45) - Scattered_light_disk2.set_inclination(70) - Scattered_light_disk2.set_density_distribution({'name':'2PowerLaws', - 'ain':2,'aout':-5,'a':20, - 'e':0,'ksi0':1., - 'gamma':2.,'beta':1.}) - Scattered_light_disk2.set_phase_function({'name':'DoubleHG', - 'g':[0.6,-0.6],'weight':0.7, - 'polar':False}) - Scattered_light_disk2.print_info() - disk2 = Scattered_light_disk2.compute_scattered_light() - - # Let's turn the polarisation on now: - Scattered_light_disk2.set_phase_function({'name':'DoubleHG', - 'g':[0.6,-0.6],'weight':0.7, - 'polar':True}) - disk2_polar = Scattered_light_disk2.compute_scattered_light() - diff --git a/vip_hci/metrics/stim.py b/vip_hci/metrics/stim.py index 82a3d55c..e85c693b 100644 --- a/vip_hci/metrics/stim.py +++ b/vip_hci/metrics/stim.py @@ -10,28 +10,29 @@ | `doi:10.1093/mnras/stz1350 `_ """ __author__ = 'Benoit Pairet' -__all__ = ['compute_stim_map', - 'compute_inverse_stim_map'] +__all__ = ['stim_map', + 'inverse_stim_map'] import numpy as np from ..preproc import cube_derotate from ..var import get_circle -def compute_stim_map(cube_der): +def stim_map(cube_der): """ Computes the STIM detection map. Parameters ---------- cube_der : 3d numpy ndarray Input de-rotated cube, e.g. ``residuals_cube_`` output from - ``vip_hci.pca.pca``. + ``vip_hci.psfsub.pca``. Returns ------- detection_map : 2d ndarray STIM detection map. """ + t, n, _ = cube_der.shape mu = np.mean(cube_der, axis=0) sigma = np.sqrt(np.var(cube_der, axis=0)) @@ -40,14 +41,14 @@ def compute_stim_map(cube_der): return get_circle(detection_map, int(np.round(n/2.))) -def compute_inverse_stim_map(cube, angle_list, **rot_options): +def inverse_stim_map(cube, angle_list, **rot_options): """ Computes the inverse STIM detection map. Parameters ---------- cube : 3d numpy ndarray - Non de-rotated residuals from reduction algorithm, eg. output residuals - from ``vip_hci.pca.pca``. + Non de-rotated residuals from reduction algorithm, eg. + ``residuals_cube`` output from ``vip_hci.psfsub.pca``. angle_list : numpy ndarray, 1d Corresponding parallactic angle for each frame. rot_options: dictionary, optional @@ -58,10 +59,11 @@ def compute_inverse_stim_map(cube, angle_list, **rot_options): Returns ------- - inverse_stim_map : 2d ndarray + inv_stim_map : 2d ndarray Inverse STIM detection map. """ + t, n, _ = cube.shape cube_inv_der = cube_derotate(cube, -angle_list, **rot_options) - inverse_stim_map = compute_stim_map(cube_inv_der) - return inverse_stim_map \ No newline at end of file + inv_stim_map = stim_map(cube_inv_der) + return inv_stim_map \ No newline at end of file diff --git a/vip_hci/negfc/__init__.py b/vip_hci/negfc/__init__.py deleted file mode 100644 index 2bc45b43..00000000 --- a/vip_hci/negfc/__init__.py +++ /dev/null @@ -1,30 +0,0 @@ -r""" -Subpackge ``negfc`` contains an ensemble of algorithms for planet position and -flux estimation. It works with as an optimization procedure using a Simplex -approach or coupled with Monte Carlo methods for posterior sampling. Thanks to -this we can get the posterior distributions of the positions and fluxes allowing -us to define proper error-bars for the photometry and position of the companion. -2 possible ways of sampling the posteriors are possible: using ``emcee`` and its -Affine Invariant MCMC or ``nestle`` with a single ellipsoid nested sampling -procedure (much faster). - -The main idea of the NegFC is to inject negative fake companions (candidates) -with varying position and flux in order to minimize (in the case of the Simplex) -a function of merit. This function of merit is defined as: - -$chi^2 = sum(\|I_j\|)$, - -where $j \in {1,...,N}$ and $N$ the total number of pixels contained in a -circular aperture (4xfwhm) centered on the injection position. This $chi^2$ is -measured on the PCA-processed frame or cube of residuals. - -""" - - -from .simplex_fmerit import * -from .mcmc_sampling import * -from .nested_sampling import * -from .simplex_optim import * -from .speckle_noise import * -from .utils_mcmc import * -from .utils_negfc import * \ No newline at end of file diff --git a/vip_hci/preproc/subsampling.py b/vip_hci/preproc/subsampling.py index 240bfbba..07144f5a 100644 --- a/vip_hci/preproc/subsampling.py +++ b/vip_hci/preproc/subsampling.py @@ -58,15 +58,13 @@ def cube_collapse(cube, mode='median', n=50, w=None): frame = np.nanmax(arr, axis=0) elif mode == 'trimmean': N = arr.shape[0] - if N % 2 == 0: - k = (N - n)//2 - else: - k = (N - n)/2 - + k = (N - n)//2 + if N%2 != n%2: + n+=1 frame = np.empty_like(arr[0]) for index, _ in np.ndenumerate(arr[0]): sort = np.sort(arr[:, index[0], index[1]]) - frame[index] = np.nanmean(sort[k:N-k]) + frame[index] = np.nanmean(sort[k:k+n]) elif mode == 'wmean': arr[np.where(np.isnan(arr))]=0 # to avoid product with nan frame = np.inner(w, np.moveaxis(arr,0,-1)) diff --git a/vip_hci/var/fit_2d.py b/vip_hci/var/fit_2d.py index b852c7de..2c5dc1ea 100644 --- a/vip_hci/var/fit_2d.py +++ b/vip_hci/var/fit_2d.py @@ -14,10 +14,12 @@ import numpy as np import pandas as pd -try: +from packaging import version +import photutils +if version.parse(photutils.__version__) >= version.parse('0.3'): # for photutils version >= '0.3' use photutils.centroids.centroid_com from photutils.centroids import centroid_com as cen_com -except: +else: # for photutils version < '0.3' use photutils.centroid_com import photutils.centroid_com as cen_com from hciplot import plot_frames diff --git a/vip_hci/var/shapes.py b/vip_hci/var/shapes.py index 8dba3387..2c67cb29 100644 --- a/vip_hci/var/shapes.py +++ b/vip_hci/var/shapes.py @@ -274,8 +274,8 @@ def get_circle(array, radius, cy=None, cx=None, mode="mask"): Notes ----- An alternative implementation would use ``skimage.draw.disk``. ``disk`` - performs better on large ``array``s (e.g. 1000px, 10.000px), while the - current implementation is faster for small ``array``s (e.g. 100px). See + performs better on large arrays (e.g. 1000px, 10.000px), while the + current implementation is faster for small arrays (e.g. 100px). See `test_shapes.py` for benchmark details. """