Commit 6eb717f6 authored by rfetick's avatar rfetick
Browse files

minor update on energy computation, set v1.3.2

parent 2a316505
......@@ -13,8 +13,8 @@ if _sys.version_info[0]<3:
_warnings.warn("MAOPPY was developped on Python 3, but your Python version is anterior. We hope everything will be fine")
# MAOPPY release version
__version__ = "1.3.1"
__version__ = "1.3.2"
__author__ = "Romain JL. Fétick (LAM, France)"
__date__ = "July 11th 2020"
__date__ = "May 30th 2021"
from . import utils, instrument, psfmodel
# -*- coding: utf-8 -*-
"""
Created on Sun May 30 12:21:19 2021
There is a difference between the theoretical and numerical values
of the sigma² parameter. This is due to numerical sampling and
the forced central frequency to be null.
@author: rfetick
"""
import matplotlib.pyplot as plt
import numpy as np
from maoppy.psfmodel import Psfao
from maoppy.instrument import muse_nfm
npix = 128 # pixel size of PSF
wvl = 600*1e-9 # wavelength [m]
#%% Initialize PSF model
samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((npix,npix),system=muse_nfm,samp=samp)
#%% Generate a MUSE-NFM PSF
r0 = 0.15 # Fried parameter [m]
b = 0 # Phase PSD background [rad² m²]
amp = 1.0 # Phase PSD Moffat amplitude [rad²]
alpha = np.logspace(-3,0,num=50) # Phase PSD Moffat alpha [1/m]
ratio = 1.0
theta = 0.0
beta = 1.6 # Phase PSD Moffat beta power law
sig2num = np.zeros_like(alpha)
f2D = Pmodel._fxfy
F2 = f2D[0]**2. + f2D[1]**2.
Fao = Pmodel.system.Nact/(2.0*Pmodel.system.D)
df = Pmodel._pix2freq
mask = F2 < Fao**2
for i in range(len(alpha)):
param = [r0,b,amp,alpha[i],ratio,theta,beta]
psd,integral = Pmodel.psd(param)
sig2num[i] = np.sum(psd*mask)*df**2
#%% PLOT
plt.figure(1)
plt.clf()
plt.semilogx(alpha/df,sig2num,label='numerical sum')
plt.semilogx(alpha/df,alpha*0 + amp, label='required')
plt.grid(which='both')
plt.xlabel("$\\alpha$/df")
plt.ylabel("$\\sigma^2$ [rad²]")
plt.legend()
\ No newline at end of file
......@@ -691,9 +691,10 @@ class Psfao(ParametricPSFfromPSD):
ay = alpha/ratio
param = (ax,ay,theta,beta,0,0)
moff = moffat(f2D,param,norm=Fao,removeInside=self._pix2freq/2) * (F2 < Fao**2.)
removeInside = (1+_np.sqrt(2))/2 * self._pix2freq/2 # remove central pixel in energy computation
moff = moffat(f2D,param,norm=Fao,removeInside=removeInside) * (F2 < Fao**2.)
moff[nx0,ny0] = 0.0 # Set Moffat PSD = 0 at null frequency
moff = moff / (_np.sum(moff)*self._pix2freq**2) # normalize moffat numerically to get correct A=sigma² in the AO zone
#moff = moff / (_np.sum(moff)*self._pix2freq**2) # normalize moffat numerically to get correct A=sigma² in the AO zone
# Warning: Moffat numerical normalization generates strehlOTF jump when "_k" is changed
PSD += (F2 < Fao**2.) * _np.abs(C + A*moff) # before I wrote "moffat(f2D,param,norm=Fao)" instead of "moff", so no numerical normalization
......
......@@ -9,7 +9,7 @@ Created on Mon May 27 17:11:01 2019
from setuptools import setup, find_packages
setup(name='MAOPPY',
version='1.3.1',
version='1.3.2',
url='https://gitlab.lam.fr/lam-grd-public/maoppy.git',
license='See LICENSE file',
author='Romain JL Fetick (LAM, Marseille, France)',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment