COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
shesha.init.geom_init Namespace Reference

Initialization of the system geometry and of the Telescope object. More...

Functions

def tel_init (carmaWrap_context context, conf.Param_geom p_geom, conf.Param_tel p_tel, r0=None, ittime=None, p_wfss=None, dm=None)
 Initialize the overall geometry of the AO system, including pupil and WFS. More...
 
def init_wfs_geom (conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, verbose=1)
 Compute the geometry of WFSs: valid subaps, positions of the subaps, flux per subap, etc... More...
 
def init_wfs_size (conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, verbose=1)
 Compute all the parameters usefull for further WFS image computation (array sizes) More...
 
def compute_nphotons (wfs_type, ittime, optthroughput, diam, cobs=0, nxsub=0, zerop=0, gsmag=0, lgsreturnperwatt=0, laserpower=0, verbose=1)
 Determines the number of photons TBC. More...
 
def init_pyrhr_geom (conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, bool verbose=True)
 Compute the geometry of PYRHR WFSs: valid subaps, positions of the subaps, flux per subap, etc... More...
 
def init_sh_geom (conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, bool verbose=True)
 Compute the geometry of SH WFSs: valid subaps, positions of the subaps, flux per subap, etc... More...
 
def geom_init (conf.Param_geom p_geom, conf.Param_tel p_tel, padding=2)
 Initialize the system geometry. More...
 
def geom_init_generic (p_geom, pupdiam, t_spiders=0.01, spiders_type="six", xc=0, yc=0, real=0, cobs=0)
 Initialize the system geometry. More...
 
def pad_array (A, N)
 

Detailed Description

Initialization of the system geometry and of the Telescope object.

Author
COMPASS Team https://github.com/ANR-COMPASS
Version
5.0.0
Date
2020/05/18

This file is part of COMPASS https://anr-compass.github.io/compass/

Copyright (C) 2011-2019 COMPASS Team https://github.com/ANR-COMPASS All rights reserved. Distributed under GNU - LGPL

COMPASS is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version.

COMPASS: End-to-end AO simulation tool using GPU acceleration The COMPASS platform was designed to meet the need of high-performance for the simulation of AO systems.

The final product includes a software package for simulating all the critical subcomponents of AO, particularly in the context of the ELT and a real-time core based on several control approaches, with performances consistent with its integration into an instrument. Taking advantage of the specific hardware architecture of the GPU, the COMPASS tool allows to achieve adequate execution speeds to conduct large simulation campaigns called to the ELT.

The COMPASS platform can be used to carry a wide variety of simulations to both testspecific components of AO of the E-ELT (such as wavefront analysis device with a pyramid or elongated Laser star), and various systems configurations such as multi-conjugate AO.

COMPASS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with COMPASS. If not, see https://www.gnu.org/licenses/lgpl-3.0.txt.

Function Documentation

◆ compute_nphotons()

def shesha.init.geom_init.compute_nphotons (   wfs_type,
  ittime,
  optthroughput,
  diam,
  cobs = 0,
  nxsub = 0,
  zerop = 0,
  gsmag = 0,
  lgsreturnperwatt = 0,
  laserpower = 0,
  verbose = 1 
)

Determines the number of photons TBC.

:parameters: wfs_type (scons.WFSType) : wfs type: SH or PYRHR.

ittime (float) : 1/loop frequency [s].

optthroughput (float) : wfs global throughput.

diam (float) : telescope diameter.

cobs (float) : (optional for SH) telescope central obstruction.

nxsub (int) : (optional for PYRHR) linear number of subaps.

zerop (float) : (optional for LGS) detector zero point expressed in ph/m**2/s in the bandwidth of the WFS.

gsmag (float) : (optional for LGS) magnitude of guide star.

lgsreturnperwatt (float) : (optional for NGS) return per watt factor (high season : 10 ph/cm2/s/W).

laserpower (float) : (optional for NGS) laser power in W.

verbose (bool) : (optional) display informations if True.

for PYRHR WFS: nphotons = compute_nphotons(scons.WFSType.PYRHR, ittime, optthroughput, diam, cobs=?, zerop=?, gsmag=?) for NGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime, optthroughput, diam, nxsub=?, zerop=?, gsmag=?) for LGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime, optthroughput, diam, nxsub=?, lgsreturnperwatt=?, laserpower=?)

Definition at line 376 of file geom_init.py.

376  lgsreturnperwatt=?, laserpower=?)
377  '''
378  surface = 0
379  nphotons = 0
380  if (wfs_type == scons.WFSType.PYRHR or wfs_type == scons.WFSType.PYRLR):
381  surface = np.pi / 4. * (1 - cobs**2.) * diam**2.
382  elif (wfs_type == scons.WFSType.SH):
383  # from the guide star
384  if (laserpower == 0):
385  if (zerop == 0):
386  zerop = 1e11
387  surface = (diam / nxsub)**2.
388  # include throughput to WFS for unobstructed
389  # subaperture per iteration
390  else: # we are dealing with a LGS
391  nphotons = lgsreturnperwatt * laserpower * \
392  optthroughput * (diam / nxsub) ** 2. * 1e4 * ittime
393  # detected by WFS
394  # ... for given power include throughput to WFS
395  # for unobstructed subaperture per iteration
396  if (verbose):
397  print("nphotons : ", nphotons)
398  return nphotons
399  else:
400  raise RuntimeError("WFS unknown")
401 
402  nphotons = zerop * 10. ** (-0.4 * gsmag) * ittime * \
403  optthroughput * surface
404 
405  if (verbose):
406  print("nphotons : ", nphotons)
407  return nphotons
408 
409 
Here is the caller graph for this function:

◆ geom_init()

def shesha.init.geom_init.geom_init ( conf.Param_geom  p_geom,
conf.Param_tel  p_tel,
  padding = 2 
)

Initialize the system geometry.

:parameters: p_geom (Param_geom) : geometry settings p_tel (Param_tel) : telescope settings padding (optional) : padding factor for PYRHR geometry

Definition at line 821 of file geom_init.py.

821  """
822 
823  # First power of 2 greater than pupdiam
824  p_geom.ssize = int(2**np.ceil(np.log2(p_geom.pupdiam) + 1))
825  # Using images centered on 1/2 pixels
826  p_geom.cent = p_geom.ssize / 2 - 0.5
827 
828  p_geom._p1 = int(np.ceil(p_geom.cent - p_geom.pupdiam / 2.))
829  p_geom._p2 = int(np.floor(p_geom.cent + p_geom.pupdiam / 2.))
830 
831  p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
832 
833  p_geom._n = p_geom.pupdiam + 2 * padding
834  p_geom._n1 = p_geom._p1 - padding
835  p_geom._n2 = p_geom._p2 + padding
836 
837  cent = p_geom.pupdiam / 2. - 0.5
838 
839  # Useful pupil
840  p_geom._spupil = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
841  cent).astype(np.float32)
842 
843  # large pupil (used for image formation)
844  p_geom._ipupil = util.pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
845 
846  # useful pupil + 4 pixels
847  p_geom._mpupil = util.pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
848 
849  if (p_tel.std_piston and p_tel.std_tt):
850  p_geom._phase_ab_M1 = mkP.make_phase_ab(p_geom.pupdiam, p_geom.pupdiam, p_tel,
851  p_geom._spupil, cent,
852  cent).astype(np.float32)
853  p_geom._phase_ab_M1_m = util.pad_array(p_geom._phase_ab_M1,
854  p_geom._n).astype(np.float32)
855  #TODO: apodizer
856  """
857  if p_geom.apod:
858  if p_geom.apodFile is None or p_geom.apodFile == '':
859  apod_filename = shesha_savepath + \
860  "apodizer/SP_HARMONI_I4_C6_N1024.npy"
861  p_geom._apodizer = makeA.make_apodizer(
862  p_geom.pupdiam, p_geom.pupdiam,
863  apod_filename.encode(), 180. / 12.).astype(np.float32)
864  else:
865  """
866  p_geom._apodizer = np.ones(p_geom._spupil.shape, dtype=np.int32)
867  p_geom._pixsize = p_tel.diam / p_geom.pupdiam
868  p_geom.is_init = True
869 
870 
Here is the caller graph for this function:

◆ geom_init_generic()

def shesha.init.geom_init.geom_init_generic (   p_geom,
  pupdiam,
  t_spiders = 0.01,
  spiders_type = "six",
  xc = 0,
  yc = 0,
  real = 0,
  cobs = 0 
)

Initialize the system geometry.

:parameters: pupdiam (long) : linear size of total pupil

t_spiders (float) : secondary supports ratio.

spiders_type (str) : secondary supports type: "four" or "six".

xc (int)

yc (int)

real (int)

cobs (float) : central obstruction ratio.

Definition at line 888 of file geom_init.py.

888  cobs: (float) : central obstruction ratio.
889  """
890  # Initialize the system pupil
891  # first poxer of 2 greater than pupdiam
892  p_geom.ssize = int(2**np.ceil(np.log2(pupdiam) + 1))
893  # using images centered on 1/2 pixels
894  p_geom.cent = p_geom.ssize / 2 + 0.5
895  # valid pupil geometry
896  pupdiam = int(pupdiam)
897  p_geom._p1 = int(np.ceil(p_geom.cent - pupdiam / 2.))
898  p_geom._p2 = int(np.floor(p_geom.cent + pupdiam / 2.))
899  p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
900  p_geom._n = p_geom.pupdiam + 4
901  p_geom._n1 = p_geom._p1 - 2
902  p_geom._n2 = p_geom._p2 + 2
903 
904  # useful pupil
905  p_geom._spupil = mkP.make_pupil_generic(pupdiam, pupdiam, t_spiders, spiders_type,
906  xc, yc, real, cobs)
907 
908  # large pupil (used for image formation)
909  p_geom._ipupil = pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
910 
911  # useful pupil + 4 pixels
912  p_geom._mpupil = pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
913 
914 
Here is the call graph for this function:

◆ init_pyrhr_geom()

def shesha.init.geom_init.init_pyrhr_geom ( conf.Param_wfs  p_wfs,
float  r0,
conf.Param_tel  p_tel,
conf.Param_geom  p_geom,
float  ittime,
bool   verbose = True 
)

Compute the geometry of PYRHR WFSs: valid subaps, positions of the subaps, flux per subap, etc...

:parameters: p_wfs (Param_wfs) : wfs settings

r0 (float) : atmos r0 @ 0.5 microns

p_tel (Param_tel) : telescope settings

geom (Param_geom) : geom settings

ittime (float) : 1/loop frequency [s]

verbose (bool) : (optional) display informations if True

Definition at line 427 of file geom_init.py.

427 
428  """
429 
430  p_wfs.npix = p_wfs._pdiam
431  p_wfs._Ntot = p_geom._n
432 
433  # Creating field stop mask
434  fsradius_pixels = int(p_wfs.fssize / p_wfs._qpixsize / 2.)
435  if (p_wfs.fstop == scons.FieldStopType.ROUND):
436  focmask = util.dist(p_wfs._Nfft, xc=p_wfs._Nfft / 2. - 0.5,
437  yc=p_wfs._Nfft / 2. - 0.5) < (fsradius_pixels)
438  elif (p_wfs.fstop == scons.FieldStopType.SQUARE):
439  X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1 # TODO: +1 ??
440  x = X[1] - (p_wfs._Nfft + 1.) / 2.
441  y = X[0] - (p_wfs._Nfft + 1.) / 2.
442  focmask = (np.abs(x) <= (fsradius_pixels)) * \
443  (np.abs(y) <= (fsradius_pixels))
444  else:
445  msg = "PYRHR wfs fstop must be FieldStopType.[ROUND|SQUARE]"
446  raise ValueError(msg)
447 
448  pyr_focmask = focmask * 1.0 # np.fft.fftshift(focmask*1.0)
449  p_wfs._submask = np.fft.fftshift(pyr_focmask)
450 
451  # Creating pyramid mask
452  pyrsize = p_wfs._Nfft
453  cobs = p_tel.cobs
454  rpup = p_geom.pupdiam / 2.0
455  dpup = p_geom.pupdiam
456  nrebin = p_wfs._nrebin
457  fracsub = p_wfs.fracsub
458  if p_wfs.pyr_pup_sep == -1:
459  pup_sep = int(p_wfs.nxsub)
460  else:
461  pup_sep = p_wfs.pyr_pup_sep
462  # number of pix between two centers two pupil images
463 
464  y = np.tile(np.arange(pyrsize) - pyrsize / 2, (pyrsize, 1))
465  x = y.T
466 
467  Pangle = pup_sep * nrebin # Pyramid angle in HR pixels
468  if p_wfs.nPupils == 0:
469  p_wfs.nPupils = 4
470  # Centers is a nPupils x 2 array describing the position of the quadrant centers
471  centers = Pangle / np.sin(
472  np.pi / p_wfs.nPupils) * np.c_[np.cos(
473  (2 * np.arange(p_wfs.nPupils) + 1) * np.pi / p_wfs.nPupils),
474  np.sin((2 * np.arange(p_wfs.nPupils) + 1) *
475  np.pi / p_wfs.nPupils)]
476  # In case nPupils == 4, we put the centers back in the normal A-B-C-D ordering scheme, for misalignment processing.
477  if p_wfs.nPupils == 4:
478  centers = np.round(centers[[1, 0, 2, 3], :]).astype(np.int32)
479  # Add in the misalignments of quadrant centers, in LR px units
480  if p_wfs.pyr_misalignments is not None:
481  mis = np.asarray(p_wfs.pyr_misalignments) * nrebin
482  else:
483  mis = np.zeros((p_wfs.nPupils, 2), dtype=np.float32)
484  # Pyramid as minimal intersect of 4 tilting planes
485  pyr = 2 * np.pi / pyrsize * np.min(
486  np.asarray([(c[0] + m[0]) * x + (c[1] + m[1]) * y
487  for c, m in zip(centers, mis)]), axis=0)
488  p_wfs._halfxy = np.fft.fftshift(pyr.T)
489 
490  # Valid pixels identification
491  # Generate buffer with pupil at center
492  pup = np.zeros((pyrsize, pyrsize), dtype=np.int32)
493  pup[pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2,
494  pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2] = \
495  p_geom._mpupil
496 
497  # Shift the mask to obtain the geometrically valid pixels
498  for qIdx in range(p_wfs.nPupils):
499  quadOnCenter = np.roll(pup, tuple((centers[qIdx] + mis[qIdx]).astype(np.int32)),
500  (0, 1))
501  mskRebin = util.rebin(quadOnCenter.copy(),
502  [pyrsize // nrebin, pyrsize // nrebin])
503  if qIdx == 0:
504  stackedSubap = np.roll(mskRebin >= fracsub,
505  tuple((-centers[qIdx] / nrebin).astype(np.int32)),
506  (0, 1))
507  mskRebTot = (mskRebin >= fracsub).astype(np.int32)
508  else:
509  stackedSubap += np.roll(mskRebin >= fracsub,
510  tuple((-centers[qIdx] / nrebin).astype(np.int32)),
511  (0, 1))
512  mskRebTot += (mskRebin >= fracsub).astype(np.int32) * (qIdx + 1)
513 
514  validRow, validCol = [], []
515  # If n == 4, we need to tweak -again- the order to feed
516  # compass XY controllers in the appropriate order
517  if p_wfs.nPupils == 4:
518  centers = centers[[2, 1, 3, 0], :]
519  # mis = mis[[2, 1, 3, 0], :] # We don't need mis anymore - but if so keep the order of centers
520  for qIdx in range(p_wfs.nPupils):
521  tmpWh = np.where(
522  np.roll(stackedSubap, tuple((centers[qIdx] / nrebin).astype(np.int32)),
523  (0, 1)))
524  validRow += [tmpWh[0].astype(np.int32)]
525  validCol += [tmpWh[1].astype(np.int32)]
526  nvalid = validRow[0].size
527  validRow = np.asarray(validRow).flatten()
528  validCol = np.asarray(validCol).flatten()
529  stack, index = np.unique(np.c_[validRow, validCol], axis=0, return_index=True)
530 
531  p_wfs._nvalid = nvalid
532  p_wfs._validsubsx = validRow[np.sort(index)]
533  p_wfs._validsubsy = validCol[np.sort(index)]
534  p_wfs._hrmap = mskRebTot.astype(np.int32)
535 
536  if (p_wfs.pyr_pos is None):
537  pixsize = (np.pi * p_wfs._qpixsize) / (3600 * 180)
538  # scale_fact = 2 * np.pi / npup * \
539  # (p_wfs.Lambda / p_tel.diam / 4.848) / pixsize * p_wfs.pyr_ampl
540  # Proposition de Flo
541  scale_fact = 2 * np.pi / pyrsize * \
542  (p_wfs.Lambda * 1e-6 / p_tel.diam) / pixsize * p_wfs.pyr_ampl
543  cx = scale_fact * \
544  np.sin((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
545  cy = scale_fact * \
546  np.cos((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
547  p_wfs.set_pyr_scale_pos(scale_fact)
548  # mod_npts = p_wfs.pyr_npts #UNUSED
549  else:
550  if (verbose):
551  print("Using user-defined positions [arcsec] for the pyramid modulation")
552  cx = p_wfs.pyr_pos[:, 0] / p_wfs._qpixsize
553  cy = p_wfs.pyr_pos[:, 1] / p_wfs._qpixsize
554  # mod_npts=cx.shape[0] #UNUSED
555 
556  p_wfs._pyr_cx = cx.copy()
557  p_wfs._pyr_cy = cy.copy()
558 
559  # telSurf = np.pi / 4. * (1 - p_tel.cobs**2.) * p_tel.diam**2.
560  # p_wfs._nphotons = p_wfs.zerop * \
561  # 10. ** (-0.4 * p_wfs.gsmag) * ittime * \
562  # p_wfs.optthroughput * telSurf
563  p_wfs._nphotons = compute_nphotons(scons.WFSType.PYRHR, ittime, p_wfs.optthroughput,
564  p_tel.diam, cobs=p_tel.cobs, zerop=p_wfs.zerop,
565  gsmag=p_wfs.gsmag, verbose=verbose)
566 
567  # spatial filtering by the pixel extent:
568  # *2/2 intended. min should be 0.40 = sinc(0.5)^2.
569  y = np.tile(np.arange(pyrsize) - pyrsize // 2, (pyrsize, 1))
570  x = y.T
571  x = x * 1. / pyrsize
572  y = y * 1. / pyrsize
573  sincar = np.fft.fftshift(np.sinc(x) * np.sinc(y))
574 
575  # sincar = np.roll(np.pi*x*np.pi*y,x.shape[1],axis=1)
576  p_wfs._sincar = sincar.astype(np.float32)
577 
578  #pup = p_geom._mpupil
579  a = pyrsize // nrebin
580  b = p_geom._n // nrebin
581  pupvalid = stackedSubap[a // 2 - b // 2:a // 2 + b // 2, a // 2 - b // 2:a // 2 +
582  b // 2]
583  p_wfs._isvalid = pupvalid.T.astype(np.int32)
584 
585  validsubsx = np.where(pupvalid)[0].astype(np.int32)
586  validsubsy = np.where(pupvalid)[1].astype(np.int32)
587 
588  istart = np.arange(p_wfs.nxsub + 2) * p_wfs.npix
589  jstart = np.copy(istart)
590 
591  # sorting out valid subaps
592  fluxPerSub = np.zeros((p_wfs.nxsub + 2, p_wfs.nxsub + 2), dtype=np.float32)
593 
594  for i in range(p_wfs.nxsub + 2):
595  indi = istart[i]
596  for j in range(p_wfs.nxsub + 2):
597  indj = jstart[j]
598  fluxPerSub[i, j] = np.sum(
599  p_geom._mpupil[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix])
600 
601  fluxPerSub = fluxPerSub / p_wfs.nxsub**2.
602 
603  p_wfs._fluxPerSub = fluxPerSub.copy()
604 
605  phasemap = np.zeros((p_wfs.npix * p_wfs.npix, p_wfs._nvalid), dtype=np.int32)
606  X = np.indices((p_geom._n, p_geom._n)) # we need c-like indice
607  tmp = X[1] + X[0] * p_geom._n
608 
609  pyrtmp = np.zeros((p_geom._n, p_geom._n), dtype=np.int32)
610 
611  for i in range(len(validsubsx)):
612  indi = istart[validsubsy[i]] # +2-1 (yorick->python
613  indj = jstart[validsubsx[i]]
614  phasemap[:, i] = tmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix].flatten("C")
615  pyrtmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix] = i
616 
617  p_wfs._phasemap = phasemap
618 
619  p_wfs._pyr_offsets = pyrtmp # pshift
620 
621 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_sh_geom()

def shesha.init.geom_init.init_sh_geom ( conf.Param_wfs  p_wfs,
float  r0,
conf.Param_tel  p_tel,
conf.Param_geom  p_geom,
float  ittime,
bool   verbose = True 
)

Compute the geometry of SH WFSs: valid subaps, positions of the subaps, flux per subap, etc...

:parameters: p_wfs (Param_wfs) : wfs settings

r0 (float) : atmos r0 @ 0.5 microns

p_tel (Param_tel) : telescope settings

geom (Param_geom) : geom settings

ittime (float) : 1/loop frequency [s]

verbose (bool) : (optional) display informations if True

Definition at line 639 of file geom_init.py.

639 
640  """
641  p_wfs.nPupils = 1
642  # this is the i,j index of lower left pixel of subap in _spupil
643  istart = ((np.linspace(0, p_geom.pupdiam, p_wfs.nxsub + 1))[:-1]).astype(np.int64)
644  # Translation in _mupil useful for raytracing
645  istart += 2
646  jstart = np.copy(istart)
647 
648  # sorting out valid subaps
649  fluxPerSub = np.zeros((p_wfs.nxsub, p_wfs.nxsub), dtype=np.float32)
650 
651  for i in range(p_wfs.nxsub):
652  indi = istart[i] # +2-1 (yorick->python)
653  for j in range(p_wfs.nxsub):
654  indj = jstart[j] # +2-1 (yorick->python)
655  fluxPerSub[i, j] = np.sum(
656  p_geom._mpupil[indi:indi + p_wfs._pdiam, indj:indj + p_wfs._pdiam])
657  # fluxPerSub[i,j] = np.where(p_geom._mpupil[indi:indi+pdiam,indj:indj+pdiam] > 0)[0].size
658 
659  fluxPerSub = fluxPerSub / p_wfs._pdiam**2.
660 
661  pupvalid = (fluxPerSub >= p_wfs.fracsub) * 1
662  p_wfs._isvalid = pupvalid.astype(np.int32)
663  p_wfs._nvalid = int(np.sum(p_wfs._isvalid))
664  p_wfs._fluxPerSub = fluxPerSub.copy()
665  validx = np.where(p_wfs._isvalid.T)[1].astype(np.int32)
666  validy = np.where(p_wfs._isvalid.T)[0].astype(np.int32)
667  p_wfs._validsubsx = validx
668  p_wfs._validsubsy = validy
669  p_wfs._validpuppixx = validx * p_wfs._pdiam + 2
670  p_wfs._validpuppixy = validy * p_wfs._pdiam + 2
671 
672  # this defines how we cut the phase into subaps
673  phasemap = np.zeros((p_wfs._pdiam * p_wfs._pdiam, p_wfs._nvalid), dtype=np.int32)
674 
675  X = np.indices((p_geom._n, p_geom._n))
676  tmp = X[1] + X[0] * p_geom._n
677 
678  n = p_wfs._nvalid
679  # for i in range(p_wfs._nvalid):
680  for i in range(n):
681  indi = istart[p_wfs._validsubsy[i]] # +2-1 (yorick->python)
682  indj = jstart[p_wfs._validsubsx[i]]
683  phasemap[:, i] = tmp[indi:indi + p_wfs._pdiam, indj:indj +
684  p_wfs._pdiam].flatten()
685  p_wfs._phasemap = phasemap
686  p_wfs._validsubsx *= p_wfs.npix
687  p_wfs._validsubsy *= p_wfs.npix
688 
689  # this is a phase shift of 1/2 pix in x and y
690  halfxy = np.linspace(0, 2 * np.pi, p_wfs._Nfft + 1)[0:p_wfs._pdiam] / 2.
691  halfxy = np.tile(halfxy, (p_wfs._pdiam, 1))
692  halfxy += halfxy.T
693  # p_wfs._halfxy = <float*>(halfxy*0.).data #dont work: half*0 is temp
694  # python obj
695 
696  if (p_wfs.npix % 2 == 1 and p_wfs._nrebin % 2 == 1):
697  # p_wfs._halfxy = <float*>(halfxy*0.)
698  halfxy = np.zeros((p_wfs._pdiam, p_wfs._pdiam), dtype=np.float32)
699  p_wfs._halfxy = halfxy.astype(np.float32)
700  else:
701  p_wfs._halfxy = halfxy.astype(np.float32)
702 
703  # this defines how we create a larger fov if required
704  if (p_wfs._Ntot != p_wfs._Nfft):
705  indi = int((p_wfs._Ntot - p_wfs._Nfft) / 2.) # +1 -1 (yorick>python)
706  indj = int(indi + p_wfs._Nfft)
707  X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1 # TODO: +1 ??
708  x = X[1]
709  y = X[0]
710  # hrpix
711  tmp = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
712  tmp[indi:indj, indi:indj] = np.roll(x + (y - 1) * p_wfs._Nfft, p_wfs._Nfft // 2,
713  axis=0)
714  tmp[indi:indj, indi:indj] = np.roll(tmp[indi:indj, indi:indj], p_wfs._Nfft // 2,
715  axis=1)
716  # hrmap=roll(hrpix)
717  tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=0)
718  tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=1)
719 
720  tmp = np.where(tmp.flatten())[0]
721 
722  p_wfs._hrmap = np.copy(tmp.astype(np.int32))
723  p_wfs._hrmap = np.reshape(p_wfs._hrmap, (p_wfs._hrmap.shape[0], 1))
724  # must be set even if unused
725 
726  else:
727  tmp = np.zeros((2, 2))
728  p_wfs._hrmap = np.copy(tmp.astype(np.int32))
729  # must be set even if unused
730 
731  if (p_wfs._nrebin * p_wfs.npix % 2 != p_wfs._Ntot % 2):
732  # +2-1 yorick>python
733  indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 1
734  else:
735  indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 0
736  indj = int(indi + p_wfs._nrebin * p_wfs.npix)
737 
738  X = np.indices((p_wfs._nrebin * p_wfs.npix, p_wfs._nrebin * p_wfs.npix))
739  x = (X[1] / p_wfs._nrebin).astype(np.int64)
740  y = (X[0] / p_wfs._nrebin).astype(np.int64)
741 
742  # binindices
743  binindices = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
744  binindices[indi:indj, indi:indj] = x + y * p_wfs.npix + 1
745 
746  binmap = np.zeros((p_wfs._nrebin * p_wfs._nrebin, p_wfs.npix * p_wfs.npix))
747 
748  X = np.indices((p_wfs._Ntot, p_wfs._Ntot))
749  tmp = X[1] + X[0] * p_wfs._Ntot
750 
751  if (p_wfs.gsalt <= 0):
752  binindices = np.roll(binindices, binindices.shape[0] // 2, axis=0)
753  binindices = np.roll(binindices, binindices.shape[1] // 2, axis=1)
754 
755  for i in range(p_wfs.npix * p_wfs.npix):
756  binmap[:, i] = tmp[np.where(binindices == i + 1)]
757  # binmap=np.reshape(binmap.flatten("F"),(binmap.shape[0],binmap.shape[1]))
758  p_wfs._binmap = np.copy(binmap.astype(np.int32))
759 
760  dr0 = p_tel.diam / r0 * \
761  (0.5 / p_wfs.Lambda) ** 1.2 / \
762  np.cos(p_geom.zenithangle * CONST.DEG2RAD) ** 0.6
763  fwhmseeing = p_wfs.Lambda / \
764  (p_tel.diam / np.sqrt(p_wfs.nxsub ** 2. + (dr0 / 1.5) ** 2.)) / 4.848
765  kernelfwhm = np.sqrt(fwhmseeing**2. + p_wfs.kernel**2.)
766 
767  tmp = util.makegaussian(p_wfs._Ntot, kernelfwhm / p_wfs._qpixsize, p_wfs._Ntot // 2,
768  p_wfs._Ntot // 2).astype(np.float32)
769 
770  tmp = np.roll(tmp, tmp.shape[0] // 2, axis=0)
771  tmp = np.roll(tmp, tmp.shape[1] // 2, axis=1)
772 
773  tmp[0, 0] = 1. # this insures that even with fwhm=0, the kernel is a dirac
774  tmp = tmp / np.sum(tmp)
775  tmp = np.fft.fft2(tmp).astype(np.complex64) / (p_wfs._Ntot * p_wfs._Ntot)
776  p_wfs._ftkernel = np.copy(tmp).astype(np.complex64)
777 
778  # dealing with photometry
779  # telSurf = np.pi / 4. * (1 - p_tel.cobs**2.) * p_tel.diam**2.
780 
781  # from the guide star
782  if (p_wfs.gsalt == 0):
783  if (p_wfs.zerop == 0):
784  p_wfs.zerop = 1e11
785  # p_wfs._nphotons = p_wfs.zerop * 10 ** (-0.4 * p_wfs.gsmag) * \
786  # p_wfs.optthroughput * \
787  # (p_tel.diam / p_wfs.nxsub) ** 2. * ittime
788  # include throughput to WFS
789  # for unobstructed subaperture
790  # per iteration
791  p_wfs._nphotons = compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
792  p_tel.diam, nxsub=p_wfs.nxsub,
793  zerop=p_wfs.zerop, gsmag=p_wfs.gsmag,
794  verbose=verbose)
795 
796  else: # we are dealing with a LGS
797  # p_wfs._nphotons = p_wfs.lgsreturnperwatt * \
798  # p_wfs.laserpower * \
799  # p_wfs.optthroughput * \
800  # (p_tel.diam / p_wfs.nxsub) ** 2. * 1e4 * \
801  # ittime
802  # detected by WFS
803  # ... for given power
804  # include throughput to WFS
805  # for unobstructed subaperture
806  # per iteration
807  p_wfs._nphotons = compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
808  p_tel.diam, nxsub=p_wfs.nxsub,
809  lgsreturnperwatt=p_wfs.lgsreturnperwatt,
810  laserpower=p_wfs.laserpower, verbose=verbose)
811 
812 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_wfs_geom()

def shesha.init.geom_init.init_wfs_geom ( conf.Param_wfs  p_wfs,
float  r0,
conf.Param_tel  p_tel,
conf.Param_geom  p_geom,
float  ittime,
  verbose = 1 
)

Compute the geometry of WFSs: valid subaps, positions of the subaps, flux per subap, etc...

:parameters: p_wfs (Param_wfs) : wfs settings

r0 (float) : atmos r0 @ 0.5 microns

p_tel (Param_tel) : telescope settings

geom (Param_geom) : geom settings

ittime (float) : 1/loop frequency [s]

verbose (int) : (optional) display informations if 0

Definition at line 129 of file geom_init.py.

129 
130  """
131 
132  if p_geom.pupdiam:
133  if p_wfs.type == scons.WFSType.SH or p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR:
134  pdiam = p_geom.pupdiam // p_wfs.nxsub
135  if (p_geom.pupdiam % p_wfs.nxsub > 0):
136  pdiam += 1
137 
138  else:
139  pdiam = -1
140 
141  p_wfs._pdiam = pdiam
142 
143  init_wfs_size(p_wfs, r0, p_tel, verbose)
144 
145  if not p_geom.is_init:
146  # this is the wfs with largest # of subaps
147  # the overall geometry is deduced from it
148  #if not p_geom.pupdiam:
149  if p_geom.pupdiam != 0 and p_geom.pupdiam != p_wfs._pdiam * p_wfs.nxsub:
150  print("WARNING: pupdiam set value not correct")
151  p_geom.pupdiam = p_wfs._pdiam * p_wfs.nxsub
152  print("pupdiam used: ", p_geom.pupdiam)
153  if p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR:
154  geom_init(p_geom, p_tel, padding=p_wfs._nrebin)
155  elif (p_wfs.type == scons.WFSType.SH):
156  geom_init(p_geom, p_tel)
157  else:
158  raise RuntimeError("This WFS can not be used")
159 
160  if (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
161  init_pyrhr_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
162  elif (p_wfs.type == scons.WFSType.SH):
163  init_sh_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
164  else:
165  raise RuntimeError("This WFS can not be used")
166 
167 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_wfs_size()

def shesha.init.geom_init.init_wfs_size ( conf.Param_wfs  p_wfs,
float  r0,
conf.Param_tel  p_tel,
  verbose = 1 
)

Compute all the parameters usefull for further WFS image computation (array sizes)

:parameters: p_wfs (Param_wfs) : wfs settings

r0 (float) : atmos r0 @ 0.5 microns

p_tel (Param_tel) : telescope settings

verbose (int) : (optional) display informations if 0

Scheme to determine arrays sizes sh : k = 6 p = k * d/r0 n = int(2*d*v/lambda/CONST.RAD2ARCSEC)+1 N = fft_goodsize(k*n/v*lambda/r0*CONST.RAD2ARCSEC) u = k * lambda / r0 * CONST.RAD2ARCSEC / N n = v/u - int(v/u) > 0.5 ? int(v/u)+1 : int(v/u) v = n * u Nt = v * Npix

pyr : Fs = field stop radius in arcsec N size of big array for FFT in pixels P pupil diameter in pixels D diameter of telescope in m Nssp number of pyr measurement points in the pupil

Rf = Fs . N . D / lambda / P ideally we choose : Fs = lambda / D . Nssp / 2

if we want good sampling of r0 (avoid aliasing of speckles) P > D / r0 . m with m = 2 or 3

to get reasonable space between pupil images : N > P.(2 + 3S) with S close to 1 N must be a power of 2

to ease computation lets assume that Nssp is a divider of P scaling factor between high res pupil images in pyramid model and actual pupil size on camera images would be P / Nssp

Definition at line 213 of file geom_init.py.

213  """
214 
215  r0 = r0 * (p_wfs.Lambda * 2)**(6. / 5)
216 
217  if (r0 != 0):
218  if (verbose):
219  print("r0 for WFS :", "%3.2f" % r0, " m")
220  # seeing = CONST.RAD2ARCSEC * (p_wfs.lambda * 1.e-6) / r0
221  if (verbose):
222  print("seeing for WFS : ",
223  "%3.2f" % (CONST.RAD2ARCSEC * (p_wfs.Lambda * 1.e-6) / r0), "\"")
224 
225  if (p_wfs._pdiam <= 0):
226  # this case is usualy for the wfs with max # of subaps
227  # we look for the best compromise between pixsize and fov
228  if (p_wfs.type == scons.WFSType.SH):
229 
230  subapdiam = p_tel.diam / float(p_wfs.nxsub) # diam of subap
231  k = 6
232  pdiam = int(k * subapdiam / r0) # number of phase points per subap
233  if (pdiam < 16):
234  pdiam = 16
235 
236  # Must be even to keep ssp and actuators grids aligned in the pupil
237  if ((pdiam * p_wfs.nxsub) % 2):
238  pdiam += 1
239 
240  nrebin = int(2 * subapdiam * p_wfs.pixsize /
241  (p_wfs.Lambda * 1.e-6) / CONST.RAD2ARCSEC) + 1
242  nrebin = max(2, nrebin)
243  # first atempt on a rebin factor
244 
245  # since we clipped pdiam we have to be carreful in nfft computation
246  Nfft = util.fft_goodsize(
247  int(pdiam / subapdiam * nrebin / p_wfs.pixsize * CONST.RAD2ARCSEC *
248  (p_wfs.Lambda * 1.e-6)))
249 
250  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
251  # while (pdiam % p_wfs.npix != 0) pdiam+=1;
252  k = 3 if p_wfs.type == scons.WFSType.PYRHR else 2
253  pdiam = int(p_tel.diam / r0 * k)
254  while (pdiam % p_wfs.nxsub != 0):
255  pdiam += 1 # we choose to have a multiple of p_wfs.nxsub
256  pdiam = pdiam // p_wfs.nxsub
257  if (pdiam < 8):
258  pdiam = 8
259 
260  # quantum pixel size
261  else:
262  pdiam = p_wfs._pdiam
263  # this case is for a wfs with fixed # of phase points
264  Nfft = util.fft_goodsize(2 * pdiam)
265  # size of the support in fourier domain
266 
267  # qpixsize = pdiam * \
268  # (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC / Nfft
269  # # quantum pixel size
270 
271  if (p_wfs.type == scons.WFSType.SH):
272  subapdiam = p_tel.diam / float(p_wfs.nxsub) # diam of subap
273 
274  # size of the support in fourier domain
275 
276  # qpixsize = k * (p_wfs.Lambda*1.e-6) / r0 * CONST.RAD2ARCSEC / Nfft
277  qpixsize = (pdiam * (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC) / Nfft
278 
279  # # actual rebin factor
280  if (p_wfs.pixsize / qpixsize - int(p_wfs.pixsize / qpixsize) > 0.5):
281  nrebin = int(p_wfs.pixsize / qpixsize) + 1
282  else:
283  nrebin = int(p_wfs.pixsize / qpixsize)
284 
285  # actual pixel size
286  pixsize = nrebin * qpixsize
287 
288  if (pixsize * p_wfs.npix > qpixsize * Nfft):
289  Ntot = util.fft_goodsize(int(pixsize * p_wfs.npix / qpixsize) + 1)
290  else:
291  Ntot = Nfft
292 
293  if (Ntot % 2 != Nfft % 2):
294  Ntot += 1
295 
296  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
297  pdiam = pdiam * p_wfs.nxsub
298  m = 3
299  # fft_goodsize( m * pdiam)
300  Nfft = int(2**np.ceil(np.log2(m * pdiam)))
301 
302  nrebin = pdiam // p_wfs.nxsub
303  while (Nfft % nrebin != 0):
304  nrebin += 1 # we choose to have a divisor of Nfft
305  pdiam = nrebin * p_wfs.nxsub
306  Nfft = int(2**np.ceil(np.log2(m * pdiam)))
307 
308  qpixsize = (pdiam *
309  (p_wfs.Lambda * 1.e-6) / p_tel.diam * CONST.RAD2ARCSEC) / Nfft
310 
311  Ntot = Nfft // pdiam * p_wfs.nxsub
312  pixsize = qpixsize * nrebin
313  pdiam = pdiam // p_wfs.nxsub
314 
315  p_wfs._pdiam = pdiam
316  p_wfs.pixsize = pixsize
317  p_wfs._qpixsize = qpixsize
318  p_wfs._Nfft = Nfft
319  p_wfs._Ntot = Ntot
320  p_wfs._nrebin = nrebin
321  p_wfs._subapd = p_tel.diam / p_wfs.nxsub
322 
323  if (verbose):
324  if (p_wfs.type == scons.WFSType.SH):
325  print("quantum pixsize : ", "%5.4f" % qpixsize, "\"")
326  print("simulated FoV : ", "%3.2f" % (Ntot * qpixsize), "\" x ",
327  "%3.2f" % (Ntot * qpixsize), "\"")
328  print("actual pixsize : ", "%5.4f" % pixsize)
329  print("actual FoV : ", "%3.2f" % (pixsize * p_wfs.npix), "\" x ",
330  "%3.2f" % (pixsize * p_wfs.npix), "\"")
331  print("number of phase points : ", p_wfs._pdiam)
332  print("size of fft support : ", Nfft)
333  print("size of HR spot support : ", Ntot)
334 
335  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
336  print("quantum pixsize in pyr image : ", "%5.4f" % qpixsize, "\"")
337  print("simulated FoV : ", "%3.2f" % (Nfft * qpixsize), "\" x ",
338  "%3.2f" % (Nfft * qpixsize), "\"")
339  print("number of phase points : ", p_wfs._pdiam * p_wfs.nxsub)
340  print("size of fft support : ", Nfft)
341 
342 
Here is the caller graph for this function:

◆ pad_array()

def shesha.init.geom_init.pad_array (   A,
  N 
)

Definition at line 915 of file geom_init.py.

915 def pad_array(A, N):
916  S = A.shape
917  D1 = (N - S[0]) // 2
918  D2 = (N - S[1]) // 2
919  padded = np.zeros((N, N))
920  padded[D1:D1 + S[0], D2:D2 + S[1]] = A
921  return padded
Here is the caller graph for this function:

◆ tel_init()

def shesha.init.geom_init.tel_init ( carmaWrap_context  context,
conf.Param_geom  p_geom,
conf.Param_tel  p_tel,
  r0 = None,
  ittime = None,
  p_wfss = None,
  dm = None 
)

Initialize the overall geometry of the AO system, including pupil and WFS.

:parameters: context (carmaWrap_context) : context p_geom (Param_geom) : geom settings p_tel (Param_tel) : telescope settings r0 (float) : atmos r0 @ 0.5 microns ittime (float) : 1/loop frequency [s] p_wfss (list of Param_wfs) : wfs settings dm (list of Param_dm) : (optional) dms settings [=None] :return: telescope (Telescope): Telescope object

Definition at line 64 of file geom_init.py.

64 
65  """
66  if p_wfss is not None:
67  # WFS geometry
68  nsensors = len(p_wfss)
69 
70  any_sh = [o.type for o in p_wfss].count(scons.WFSType.SH) > 0
71  # dm = None
72  if (p_wfss[0].dms_seen is None and dm is not None):
73  for i in range(nsensors):
74  if (not p_wfss[i].open_loop):
75  p_wfss[i].set_dms_seen(np.arange(len(dm), dtype=np.int32))
76 
77  # first get the wfs with max # of subaps
78  # we'll derive the geometry from the requirements in terms of sampling
79  if (any_sh):
80  indmax = np.argsort([o.nxsub for o in p_wfss
81  if o.type == scons.WFSType.SH])[-1]
82  else:
83  indmax = np.argsort([o.nxsub for o in p_wfss])[-1]
84 
85  print("*-----------------------")
86  print("Computing geometry of WFS", indmax)
87 
88  init_wfs_geom(p_wfss[indmax], r0, p_tel, p_geom, ittime, verbose=1)
89  # #do the same for other wfs
90  for i in range(nsensors):
91  if (i != indmax):
92  print("*-----------------------")
93  print("Computing geometry of WFS", i)
94  init_wfs_geom(p_wfss[i], r0, p_tel, p_geom, ittime, verbose=1)
95 
96  else:
97  geom_init(p_geom, p_tel)
98 
99  telescope = Telescope(context, p_geom._spupil.shape[0],
100  np.where(p_geom._spupil > 0)[0].size,
101  (p_geom._spupil * p_geom._apodizer).astype(np.float32),
102  p_geom._mpupil.shape[0], p_geom._mpupil)
103 
104  if (p_geom._phase_ab_M1 is not None):
105  telescope.set_phase_ab_M1(p_geom._phase_ab_M1)
106  telescope.set_phase_ab_M1_m(p_geom._phase_ab_M1_m)
107 
108  return telescope
109 
110 
Here is the call graph for this function: