anisocado.psf_utils module

Originally from file _anisocado.py

anisocado.psf_utils.airmassImpact(r0_at_zenith, zenith_distance)

<r0_at_zenith> : r0 at zenith (any unit is ok) <zenith_distance> : zenith distance (degrees)

The seeing/r0 actually perceived by a telescope depends on the zenith distance. This function converts a r0 given at zenith into the real r0 actually observed by the telescope.

anisocado.psf_utils.aliasingSpectrum(kx, ky, r0, L0, M4, dssp=0.4015)

Example:

N = 512
pixelSize = 4.2
L0 = 25.0
r0 = 0.6
wavelength = 1.65e-6
rotdegree = 10.
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, rotdegree)
W = aliasingSpectrum(kx, ky, r0, L0, M4)

anisocado.psf_utils.anisoplanaticSpectrum(Cn2h, layerAltitude, L0, offx, offy, wavelength, kx, ky, Wiener, M4)

<Cn2h>list of normalised (np.sum(Cn2h)==1.0) strengh of

turbulence of each layer (no unit)

<layerAltitude> : list of altitudes of the turbulence layers (metres) <L0> : value of the outer scale (metres) <offx> : off-axis distance of the star along X axis (arcsec) <offy> : idem along Y axis (arcsec) <wavelength> : wavelength (metres) <kx> : spatial frequency produced by

the function computeSpatialFreqArrays() (metres^-1)

<ky> : idem kx <Wiener> : turbulent spectrum from function computeWiener() <M4> : frequency domain of M4 (boolean) computed by the

function M4 = defineDmFrequencyArea(kx, ky, 0)

Computes the phase spectrum (spatial frequencies) contribution due to the anisoplanatism in the AO compensation for a source located off-axis from the SCAO guide star by some amount (offx, offy). In input, the distribution of the turbulence in altitude is given.

Result is returned with a frequency corner-centred representation.

Example:

N = 512
pixelSize = 4.2
Cn2h = [0.3, 0.2, 0.2, 0.1]
layerAltitude = [0,1000,8000,15000]
L0 = 25.0
offx, offy = (34., 12.)
wavelength = 1.65e-6
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, 0)
W = computeWiener(kx, ky, L0, 1.0)
f = anisoplanaticSpectrum(Cn2h, layerAltitude, L0, offx, offy,
                          wavelength, kx, ky, W, M4)

anisocado.psf_utils.clean_psf(psf, threshold)

anisocado.psf_utils.computeBpSpectrum(kx, ky, V, Fe, tret, gain, Wiener, M4)

Example:

N = 512
pixelSize = 4.2
L0 = 25.0
r0 = 0.6
wavelength = 1.65e-6
rotdegree = 10.
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, rotdegree)
W = computeWiener(kx, ky, L0, r0)

V = 10.
Fe = 500.
tret = 0.004
gain = 0.3
f = computeBpSpectrum(kx, ky, V, Fe, tret, gain, W, M4)

anisocado.psf_utils.computeEeltOTF(pup)

anisocado.psf_utils.computeSpatialFreqArrays(N, pixelSize, wavelength)

<N> : size of the output image (NxN) <pixelSize> : size of the pixels of the psf image (mas) <wavelength> : wavelength (metres)

The function returns a tuple of spatial frequencies (m^-1) in X and Y together with the pixel scale of these (the value will be useful later on in other procedures). Results of 2D arrays are returned with a Fourier corner-centred representation.

Standalone example:

N = 512               # output will be 512x512
pixelSize = 4.2       # 4.2 mas
wavelength = 1.65e-6  # metres
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)

anisocado.psf_utils.computeWiener(kx, ky, L0, r0)

<kx>spatial frequency produced by

the function computeSpatialFreqArrays() (metres^-1)

<ky> : idem kx <L0> : value of the outer scale (metres) <r0> : value of Fried parameter r0 (metres)

The function returns the 2D spectrum of the turbulence with an outer scale L0. It is expressed in rad^2.m^2.

Result is returned with a frequency corner-centred representation.

Standalone example:

N = 512               # output will be 512x512
pixelSize = 4.2       # 4.2 mas
wavelength = 1.65e-6  # metres.  H band.
L0 = 25.              # 25 metres outer scale
r0 = 0.6              # r0 = 60cm in H band
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
W = computeWiener(kx, ky, L0, r0)

anisocado.psf_utils.convertSpectrum2Dphi(W, uk)

<W>spatial spectrum to be converted into phase structure function

in rd^2.m^2

<uk> : # size of the ‘spatial frequency pixel’ in m^-1

Converts the 2D spectrum into a phase structure function Dphi. Uses Dphi(r) = $ $ (1-cos(2.pi.k.r)) W(k) d2k Computation of Dphi is in radians^2 at the wavelength of r0.

anisocado.psf_utils.core_generatePsf(Dphi, FTOtel)

Examples

::

N = 1024 pixelSize = 4.2

# atmospheric profile (old ESO profile before 2010) layerAltitude = [47., 140, 281, 562, 1125, 2250, 4500, 9000, 18000.]

# from ref. E-SPE-ESO-276-0206_atmosphericparameters Cn2h = [52.24, 2.6, 4.44, 11.60, 9.89, 2.95, 5.98, 4.30, 6] Cn2h = np.array(Cn2h) Cn2h /= np.sum(Cn2h)

# outer scale L0 = 25.0 offx, offy = (15., 20.) wavelength = 1.65e-6 rotdegree = 10.0 kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength) M4 = defineDmFrequencyArea(kx, ky, rotdegree) r0 = 0.6

# This is the turbulent spectrum …. W = computeWiener(kx, ky, L0, r0)

# And here are some of the PSF-destroyers Waniso = anisoplanaticSpectrum(Cn2h, layerAltitude, L0, offx, offy,

wavelength, kx, ky, W, M4)

Wfit = fittingSpectrum(W, M4) nmRms = 100. Wother = otherSpectrum(nmRms, M4, uk, wavelength) Dphi = convertSpectrum2Dphi(Waniso + Wfit + Wother, uk)

# Here, i need to generate a kind of PSF or telescope OTF deadSegments = 3 pup = fake_generatePupil(N, deadSegments, rotdegree, pixelSize,

wavelength)

FTOtel = computeEeltOTF(pup)

psf = core_generatePsf(Dphi, FTOtel)

print(psf.max()) window = 100 plt.imshow( psf[N//2-window:N//2+window, N//2-window:N//2+window]**0.3 )

anisocado.psf_utils.createAdHocScaoPsf(N, pixelSize, wavelengthIR, rotdegree, r0Vis, nmRms)

<N>size of the output image, that is made to match the size

of the (square) psf image to be processed. In other words, N = psf.shape[0]

<pixelSize> : size of the pixels of the psf image (mas) <wavelengthIR> : IR wavelength (metres) <rotdegree> : pupil rotation (degrees) <r0Vis> : value of r0 in the visible (metres) <nmRms> : number of nm rms affecting the psf (nm)

Do not use that function. It’s just there to create a quicky random shitty psf. It’s there because i needed a random shitty psf. So i wrote it. And it’s still there.

Example:

N = 512
pixelSise = 4.2         # mas
wavelengthIR = 1.65e-6  # metres
rotdegree = 10.
r0Vis = 0.12
nmRms = 150.
psf, pup = createAdHocScaoPsf(N, pixelSize, wavelengthIR, rotdegree,
                              r0Vis, nmRms)

anisocado.psf_utils.defineDmFrequencyArea(kx, ky, rotdegree, dactu=0.5403)

<kx>spatial frequency produced by the

function computeSpatialFreqArrays() (metres^-1)

<ky> : idem kx <rotdegree> : rotation of M4 (degrees) <dactu> : value of actuator pitch of M4 (metres)

The function returns the 2D domain of spatial frequencies which can be impacted by M4, i.e. the M4 compensation domain. Underlying assumtions are than M4 is of an hexagonal pattern, with an inter-actuator distance <dactu> set by default to 54.03 cm.

Result is returned with a frequency corner-centred representation.

Standalone example:

N = 512               # output will be 512x512
pixelSize = 4.2       # 4.2 mas
wavelength = 1.65e-6  # metres
kx, ky = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, 0)
plt.imshow(np.fft.fftshift(M4).T, origin='l')

anisocado.psf_utils.fake_generatePupil(N, deadSegments, rotdegree, pixelSize, wavelength, rng=Generator(PCG64) at 0x7F7C6E9A3120)

<N>size of the output image, that is made to match the size

of the (square) psf image to be processed. In other words, N = psf.shape[0]

<deadSegments> : number of hexa segments of M1 that are missing <rotdegree> : pupil rotation in degrees <pixelSize> : size of the pixels of the psf image (mas) <wavelength> : wavelength (metres) <rng> : optional random number generator for reproducible results

Examples:

N = 512
deadSegments = 3
rotdegree = 14.
pixelSize = 4.2
wavelength = 1.65e-6
pup = fake_generatePupil(N, deadSegments, rotdegree, pixelSize,
                         wavelength)

anisocado.psf_utils.fittingSpectrum(Wiener, M4)

<Wiener> : turbulent spectrum from function computeWiener() <M4> : frequency domain of M4 (boolean) computed by the

function M4 = defineDmFrequencyArea(kx, ky, 0)

Returns the spatial spectrum of the (so-called) fitting error, i.e. the residual phase after a full, perfect, ideal, instantaneous, super-duper, hyper-clean, theoretical compensation of M4. It is expressed in rad^2.m^2. Result is returned with a frequency corner-centred representation.

Example:

N = 512
pixelSize = 4.2
L0 = 25.0
offx, offy = (34., 12.)
wavelength = 1.65e-6
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, 0)
r0 = 0.6
W = computeWiener(kx, ky, L0, r0)
f = fittingSpectrum(W, M4)

anisocado.psf_utils.get_atmospheric_turbulence(myProfile='EsoMedian')

Returns the relative level of turbulence at a given height

Note: The np.sum(Cn2h) = 1.0 Turbulence profile have been taken from ESO-258292: “Relevant Atmospheric Parameters for E-ELT AO Analysis and Simulations”.

The following 3 turbulence profile are currently available:

• EsoQ1

  Best atmospheric conditions - First Quartile Armazones atmospheric turbulence profile

• EsoMedian [also officialEsoMedian]

  Median atmospheric conditions - Median Armazones atmospheric turbulence profile

• EsoQ4

  Worst atmospheric conditions - First Quartile Armazones atmospheric turbulence profile

Additional turbulence profiles include:

• oldEso

  The old ESO profile from before 2010. Cn2h paramateres are taken from ref. E-SPE-ESO-276-0206_atmosphericparameters

• gendron

  The Gendron profile. Short, fast. Saves CPU. Carbon efficient. The np.sum(Cn2h) = 1.0 is hyper-guaranteed here.

Parameters:

myProfilestr, optional

Profile name: [‘EsoQ1’, ‘EsoMedian’, ‘EsoQ4’, ‘oldEso’, ‘gendron’]

Returns:

layerAltitudelist of floats

[m] height of layer above ground

Cn2hlist of floats

Relative strength of turbulence

anisocado.psf_utils.get_profile_defaults(myProfile='EsoMedian')

Data taken from ESO-258292

Parameters:

myProfilestr

[EsoMedian, EsoQ1, EsoQ4, OldEso, Gendron]

Returns:

seeingfloat

[arcsec]

zen_distfloat

[deg] Zenith distance

wind_alphafloat

Multiplication factor for the wind profile

anisocado.psf_utils.hcor(freq, Fe, tret, G, BP, an=True)

* Fonction extraite de STYC le 9 Jan 2013 *

The option an=1 sets the integrator to “analog”. Doing this, an extra 1/2 frame delay is added compared to case of the numeric integrator an=0.

<tret> is the delay expressed as a time in seconds, between the end of the integration and the start of the command.

anisocado.psf_utils.otherSpectrum(nmRms, M4, uk, wavelength)

<nmRms> : number of nm rms <M4> : frequency domain of M4 (boolean) computed by the

function M4 = defineDmFrequencyArea(kx, ky, 0)

<uk> : # size of the ‘spatial frequency pixel’ in m^-1 <wavelength> : wavelength (metres)

Example:

N = 512
pixelSize = 4.2
L0 = 25.0
offx, offy = (34., 12.)
wavelength = 1.65e-6
kx, ky, uk = computeSpatialFreqArrays(N, pixelSize, wavelength)
M4 = defineDmFrequencyArea(kx, ky, 0)
nmRms = 250.
f = otherSpectrum(nmRms, M4, uk, wavelength)

anisocado.psf_utils.r0Converter(r0, lambda1, lambda2)

Converts a r0 defined at some wavelength lambda1, into a r0’ at lambda2. Lambda1 and 2 shall be in the SAME unit. Returned r0 will be in the SAME unit than the input one.

Example:

r0_K = r0Converter(0.12, 500, 2200)

anisocado.psf_utils.round_edges(kernel, edge_width=10)