From 63824b599ca0065291455afc7bc963f90a6e47ea Mon Sep 17 00:00:00 2001
From: Romain Fetick <rfetick@OSU.PYTHEAS>
Date: Tue, 6 Jun 2023 17:21:24 +0200
Subject: [PATCH] add variance computation into Psfao
---
maoppy/psfmodel.py | 26 +++++++++++++++++++++++---
1 file changed, 23 insertions(+), 3 deletions(-)
diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index e60ace5..608d838 100644
--- a/maoppy/psfmodel.py
+++ b/maoppy/psfmodel.py
@@ -42,7 +42,7 @@ class ParametricPSF:
self.bounds = ((-np.inf,)*nparam, (np.inf,)*nparam)
def __str__(self):
- s = "ParametricPSF\n"
+ s = self.__class__.__name__ + "\n"
s += "-------------\n"
s += "pixels : (%u,%u)\n"%self.npix
s += "nb param: %u"%self._nparam
@@ -625,6 +625,25 @@ class Psfao(ParametricPSFfromPSD):
return ["r0","bck","amp","alpha","ratio","theta","beta"]
+ def var_corr(self, parampsd):
+ """Return the numerical variance on the corrected area"""
+ psd,_ = self.psd(parampsd)
+ return np.sum(psd*self._maskin) * self.pix2freq**2
+
+
+ def var_halo(self, parampsd, fovnorm=False):
+ """Return the numerical variance on the halo (fitting error)"""
+ psd,var_tot = self.psd(parampsd, fovnorm=fovnorm)
+ var_corr = np.sum(psd*self._maskin) * self.pix2freq**2
+ return var_tot - var_corr
+
+
+ def var_total(self, parampsd, fovnorm=False):
+ """Return total variance halo + corrected area"""
+ _,var_tot = self.psd(parampsd, fovnorm=fovnorm)
+ return var_tot
+
+
def psd(self, parampsd, grad=False, fovnorm=False):
"""Compute the PSD of the electromagnetic phase
PSD is given in [rad²/f²] = [rad² m²]
@@ -654,8 +673,9 @@ class Psfao(ParametricPSFfromPSD):
r0, bck, amp, alpha, ratio, theta, beta = parampsd
- if alpha < pix:
- warnings.warn("`alpha < df`. Everything works fine but the value of `amp` might not correspond exactly to the PSD integral.")
+ if alpha < (2*pix):
+ # well, the factor 2 is really empirical...
+ warnings.warn("`alpha < 2*df`. Everything works fine but the value of `amp` might not correspond exactly to the PSD integral.")
# Von-Karman
psd = (r0**(-5./3.)) * self._vk
--
GitLab