HOME   |   GITHUB   |   INSTALL   |   MANUAL   |   NEWS

This manual is for yao v5.14.1. Last modified July 2021. Expand/Collapse all detailled comments.

A session at a glance...

A basic yao session

... consist of a few steps:

  1. Craft up a parameter file describing the system you want to model, say "mysystem.par"
  2. Start yorick/yao
  3. aoread,"mysystem.par"
  4. aoinit,some keywords...
  5. aoloop,some keywords...
  6. go

Let us now enter into some details.

The Basics

Is yorick installed ? All set up as per the installation instructions ? Then go to the examples directory of the yao distribution; this might be in different locations depending on how you installed yao. To determine where is it is, run the command find_examples_path() within yao. Type the following at the unix prompt:

poliahu:~ $ yorick -i yao.i
 Copyright (c) 2005.  The Regents of the University of California.
 All rights reserved.  Yorick 2.1.05x ready.  For help type 'help'
 Yao version 4.5.0, Last modified 2009oct23

You get the yorick and yao welcome messages and the yorick prompt. Alternatively, you can start a normal yorick session and then include yao at any time by typing:

> #include "yao.i"

You can double check everything is normal by typing:

> info,aoread
func aoread(parfile)

if you get this message, you are in business. If not, fix your yorick installation. Click here to show/hide notes on Yorick.

The first thing to do is to create phase screens (to simulate turbulence). Type

> create_phase_screens,2048,256,prefix="screen"

This will create N (N=long dimension/short dimension, 8 in that case) phase screens of dimension 2048x256 suitable for use by yao. It is advised to choose dimensions that are powers of 2. Depending on your platform and CPU, it may take some time (1mn or so), as this routine is absolutely not optimized. This is a one shot run. You will not need to do that everytime you run yao, as you can, and are encouraged, to use the same phase screens. You may need to run it once more to create larger phase screens if the need arises, but that's about it. The phase screens (screen1.fits to screen8.fits in the example above) will be created in the current working directory. Move them somewhere convenient (I have them in Y_USER/data=.yorick/data in my case). You will need to edit the yao parameter files to reflect the path and names of these phase screens if you put it somewhere different or used a different name (look for "atm.screen" in the parfile).

Next, we will try to run "sh6x6", which is a simple 6x6 Shack-Hartmann example. After you have edited the "sh6x6.par" file (in the examples directory) and modified the path and filename of the newly created phase screens, if needed, type:

> aoread,"sh6x6.par"
Yao, Version 4.5.0, 2009oct23
Checking parameters ...
No field stop defined for wfs 1. Setting to 'square'
wfs(1).fssize has not been set, will be forced to subap FoV
dm(1).coupling set to 0.200000
dm(1).iffile set to sh6x6-if1.fits
dm(2).iffile set to sh6x6-if2.fits

What aoread() does is (a) read, or rather include, the parameter file, which will fill the various structures containing the definitions of the WFS, DM, loop, etc... and (b) go through a simple check of the parameters to see if anything is missing or if there are incompatible assigments, in which case it will print out an error message (hopefully understandable). Otherwise, it prints out informational messages or warnings.

Then we need to initialize the system. aoinit() will do that for us. It will initialize all the arrays (pupil, etc), initialize the system pupil, the various WFS, DM, etc. It will then compute the interaction matrix, invert it and finally plot (if requested) a graphical system configuration. The amount of information you get during the aoinit is set by sim.verbose. The default verbose in sh6x6.par is 0, which means you get no feedback at all except for warnings and error messages. Let's set sim.verbose to 1 and run aoinit():

> sim.verbose=1
> aoinit,disp=1,clean=1
Checking parameters ...

Reading phase screen "~/.yorick/data/screen1.fits"
Reading phase screen "~/.yorick/data/screen2.fits"
Reading phase screen "~/.yorick/data/screen3.fits"
Reading phase screen "~/.yorick/data/screen4.fits"

Kernel FWHM for the iMat calibration = 0.364983
Pre-computing Kernels for the SH WFS
 WFS# |       Pixel sizes         | Subap. size | Number of pixels | #photons
      | Desired  Quantum  Actual  | Max  Actual | Desired   Actual | /sub/iter
 1      0.20000  0.03182  0.19093   2.04  1.91    10x10      10x10   70735.5
NGS#1 flux varies between 42795 and 70736 photon/subap/it

> Initializing DM influence functions
  >> Computing Influence functions for mirror # 1

Creating Influence function for actuator #1 2 3 4 5 6 7 8 9 10 11 12
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
36 37 38 39 40 41 42 43 44 45 46 47 48 49

  >> Storing influence functions in sh6x6-if1.fits... Done
  >> Computing Influence functions for mirror # 2

  >> Storing influence functions in sh6x6-if2.fits... Done

  DM #1: # of valid actuators: 45. (I got rid of 4 actuators after iMat)
  >> valid I.F. stored in sh6x6-if1.fits
Computing valid to extrap. matrix for DM#1

I did not find simulModeGains.fits or it did not have the right
   number of elements (should be 47), so I have generated
   a modal gain vector with ones in it.

 >> Preparing SVD and computing command matrice
Smallests 2 normalized eigenvalues = 0.014315    1.62232e-07
4 modes discarded in the inversion

Mirror #1, stackarray, 45 actuators, conjugated @ 0 m
Mirror #2, tiptilt, 2 actuators, conjugated @ 0 m
WFS # 1, hartmann (meth. 2), 32 subap., offaxis (+0.0",+0.0"), noise enabled
D/r0 (500nm) = 42.4; 5000 iterations

At this point you have initialized everything. The aoinit() keyword clean indicates that you want to start from scratch, and ignore any influence function or interaction/command matrix files on your disk. The disp=1 is to get some graphical feedback. You are now ready to run the loop:

> loop.niter = 1000
> aoloop,disp=10
NGS#1 flux varies between 42795 and 70736 photon/subap/it

> Starting loop with 1000 iterations
> go
 Iter#  Inst.Strehl  Long expo.Strehl  Time Left
     1  0.000        0.000             00:00:18.8
    51  0.407        0.393             00:00:15.1
   101  0.412        0.410             00:00:13.7
   901  0.647        0.504             00:00:01.4
   951  0.546        0.507             00:00:00.7
Saving results in /home/frigaut/yorick-2.1/share/yao/examples/sh6x6.res (ps,imav.fits)...
time(1-2) =  9.49 ms  (WF sensing)
time(2-3) =  0.03 ms  (Reset and measurement history handling)
time(3-4) =  0.01 ms  (cMat multiplication)
time(4-5) =  1.85 ms  (DM shape computation)
time(5-6) =  1.44 ms  (Target PSFs estimation)
time(6-7) =  1.53 ms  (Displays)
time(7-8) =  0.05 ms  (Circular buffers, end-of-loop printouts)
Finished on 00:26:30
69.040800 iterations/second in average

         lambda   XPos   YPos  FWHM[mas]  Strehl  E50d[mas]  #modes comp.
Star# 1    1.65    0.0    0.0       44.1   0.507      172.0         35.1
Field Avg  1.65                     44.1   0.507      172.0
Field rms                            0.0   0.000        0.0

You have ran 1000 iterations (loop.niter in sh6x6.par is larger so we changed it before starting the loop in the example above to keep your demo time reasonable -note if you want to make it larger than the initial value, you will need to re-run aoinit()). disp=10 in the call means "do your displays every 10 iterations". Depending on your graphic card, displays can be pretty expensive (time), thus it helps to display less frequently.

At any time, while the loop is running, you have access to the yorick prompt. You can type regular commands, but the most useful are:

at the completion of the requested number of iterations, go() will automatically call the function after_loop(), which outputs a number of things as shown above (starting from "Saving results..."): Some statistics on execution time and number of iterations per seconds, and the Strehl, FWHM, etc, on every "target" for which positions were specified in the parameter file. You can call after_loop() by hand at any time also, and you can call it multiple times.

The resulting average images have been saved on disk ("sh6x6-imav.fits"), together with a small postscript file ("sh6x6.ps") that contains some graphics. Strehl, FWHM and 50 percent encircled energy are available as extern variables under the name strehl, fwhm and e50 (the averaged PSFs are also available, together with the history of DM commands, DM errors and WFS measurements if the keyword savecb= has been set):

> strehl

These variables contains the values for all the images (here there is only one, but there can be an arbitrary number of positions and wavelengths at which you want to estimate the performance). This can be useful in script, as detailled below.

Finally, note that you can have access to some documentation on each function by typing help,function_name, e.g:

> help, aoinit
 /* DOCUMENT aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
    Second function of the ao serie.
    Initialize everything in preparation for the loop (aoloop).
    Follows a call to aoread, parfile.

    disp         = set to display stuff
    clean        = if set, aoinit will start fresh. *Nothing* is kept from
Tip: sometimes the document section of a function is not up to date. You can get a peek on the actual function API by using info:
> info, aoinit
func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)

The Parameter File

Here is an example of a parameter file (parfile). Comments below.

// YAO parameter file
sim.name           = "SH6x6 w/ TT mirror and WFS, full diffraction WFS";
sim.pupildiam      = 120;
sim.debug          = 0;
sim.verbose        = 0;

atm.dr0at05mic     = 42.44;  // this is r0=0.166 at 550 nm
atm.screen	   = &(Y_USER+"data/screen"+["1","2","3","4"]+".fits");
atm.layerfrac      = &([0.4,0.2,0.3,0.1]);
atm.layerspeed     = &([11.,20,29,35]);
atm.layeralt	   = &([0.,400,6000,9000]);
atm.winddir        = &([0,0,0,0]);

nwfs               = 1; // number of WFSs (>1 if e.g. mcao)
wfs                = array(wfss,nwfs);

wfs(1).type	   = "hartmann";
wfs(1).lambda	   = 0.65;
wfs(1).gspos       = [0.,0.];
wfs(1).gsalt       = 0.;
wfs(1).gsmag       = 5.;
wfs(1).shmethod	   = 2;
wfs(1).shnxsub	   = 6;
wfs(1).pixsize     = 0.2;
wfs(1).npixels     = 10;
wfs(1).noise       = 1;
wfs(1).ron         = 3.5;
wfs(1).shthreshold = 0.;
wfs(1).nintegcycles= 1;
wfs(n).svipc       = 2;

ndm                = 2;
dm                 = array(dms,ndm);

n  =1;
dm(n).type	   = "stackarray";
dm(n).nxact	   = 7;
dm(n).pitch	   = 20;
dm(n).alt          = 0.;
dm(n).unitpervolt  = 0.01;
dm(n).push4imat    = 100;

n  =2;
dm(n).type	   = "tiptilt";
dm(n).alt          = 0.;
dm(n).unitpervolt  = 0.0005;
dm(n).push4imat    = 400;

mat.condition      = &([15.]);

tel.diam	   = 7.9;
tel.cobs	   = 0.1125;

target.lambda	   = &([1.65]);
target.xposition   = &([0.]);
target.yposition   = &([0]);
target.dispzoom    = &([1.]);

gs.zeropoint	   = 1e11;

loop.gain	   = 0.6;
loop.framedelay    = 1;
loop.niter	   = 5000;
loop.ittime	   = 2e-3;
loop.startskip     = 10;
loop.skipevery     = 10000;
loop.skipby        = 10000;
loop.modalgainfile = "simulModeGains.fits";


The parfile defines entirely your system. It consists of nine subsections in which the yao structures are filled. All of the structure members are listed in the yao structure section, together with a short explanation of each parameters (structure member). How to combine these parameters to get yao to do what you want is explained in the controlling features section below.

Several comments:

Controlling Features

Overall Geometry

In yao, you have to define the AO system you want to simulate. It starts by defining an entrance aperture (the system pupil). This is done through 2 parameters: the physical pupil size (e.g. diameter) in real world units, meters. And because yao is a monte-carlo code, that uses arrays to generate phases and PSFs, we will need a pupil array and thus a pupil diameter in pixel.

sim.pupildiam  <- this is the pupil diameter in pixel (unitless)
tel.diam       <- this is the telescope diameter in meters

Typical geometry of a system (WFS and DM)

yao also uses arrays to store the deformable mirror influence functions, etc. Generally speaking, there are two types of variable and arrays: the one referring to quantities in the near field (pupil plane or close to it, e.g. altitude layers) and the one referring to quantities in the far field (Shack-Hartmann WFS spots, PSF image, observed object, etc...).

The figure on the right illustrates how you can set up a Shack-Hartmann + Stackarray deformable mirror (allegedly the most common type of AO system to date). I will review all these parameters in sections below, but what I want to emphasize here is the following:

Field of view considerations

The size of the pupil is actually quite important. By property of the Fourier transform, the field of view in the far field will be defined by the sampling in the near field, following the relation:

FoV = lambda/ps
where ps = tel.diam/sim.pupildiam is the pixel size in meter in the near field. So, for instance,

tel.diam [m]sim.pupildiam [pixels]FoV @650nm ['']FoV @1.65microns ['']Comment
8128 2.145.44
302561.142.90Probably too low

Choose the field of view according to your application. But because eventually you will want to run simulations in presence of turbulence, I generally advise to take at least 2 points per r0 (near field), and preferably 3. For r0 = 16 cm, that means a sampling of 8 cm (resp 5.33 cm), and for a 8-m telescope that would translate into 100 points (resp 150) across the pupil (sim.pupildiam = 100). If this is not done, you will likely end up with aliasing in the SHWFS subaperture images (probably less so in the final image if it is estimated in the Near Infrared and the WFS works in the visible). This is bad and will bias your results. You have been warned.

That said, yao 5.0 has introduced a new feature for Shack-Hartmann WFS: extended field of view. When using this feature (control it with wfs.extfield), the individual subaperture imagelets are embedded into larger intermediate images before being binned and embedded into the final SHWFS spot image. This allows extending the field of view of subapertures without the prohibite cost of using a very large number of pixels in your pupil diameter (as outlined above). It is of particular interest for large LGS elongations, e.g. for ELTs. The whole process, describing the steps used to compute the SHWFS final spot image, is somewhat detailled in this page.


Historically, yao was developed assuming circular pupils. So in any case, as said above, you will have to define sim.pupildiam and tel.diam. If another pupil shape is desired, it can be defined through a call to a user defined function user_pupil(). Click here to show/hide an example.

Wavefront Sensors

Shack-Hartmann WFS

Set wfs.type="hartmann". The only other mandatory parameters to define are wfs.shnxsub, wfs.npixpersub, wfs.npixels, wfs.pixsize and wfs.lambda. Below I expand on these parameters.

Near Field

The near field parameters that need to be defined are:

For a "normal" system, we want wfs.shnxsub * wfs.npixpersub = sim.pupildiam. But as said above, you can make npixpersub larger or smaller if you want to investigate the effect of subapertures larger or smaller than ideal. Also, by setting wfs.pupoffset [meter], you can investigate the effect of misregistration of the WFS w.r.t the entrance pupil (Show/Hide details).

Far Field

The far field parameters to define for a SHWFS are:

So that of course, the subaperture field of view will be wfs.npixels * wfs.pixsize. The whole discussion about sampling above of course applies here too, and field of view can not be arbitrary large. The pixel size is also constrained by FFT properties, and yao will round the pixel size you have selected to the nearest possible value (this should be printed out on screen during aoinit if sim.verbose>=1).

Field Stops

Since v4.5.0, yao can handle field stops (before also, but not in the same easy fashion), defined using the parameters wfs.fstop, wfs.fssize and wfs.fsoffset. You have the choice between no field stop (wfs.fstop="none"), or a square or round field stop (wfs.fstop="square" or "round"). If needed, you can use your own defined geometry by setting wfs.fsname to the name of a fits file that contains the image of a field stop (show/hide how). The default is to use a square field stop of size equal to the subaperture field of view (i.e. optimal).

Subaperture Coupling

Shack-Hartmann focal plane features

Since v4.5.0, yao correctly includes coupling between neighbor subapertures: In a real system, if there is no field stop, the spot image of a subaperture can wander into a neighbor subaperture and wreak havoc the slope estimation. yao now correctly handle this behavior: all subaperture images are computed, then added to an image that will include all spots. When all spots have been computed, each spot image is extracted and the spot position is estimated from there (center of gravity, soon to come weighted CoG, etc...).

Graphic configuration

When the correct debug setting is selected (sim.debug>=1), aoinit() will display a SHWFS far field graphic configuration diagram as shown on the right (click on it to enlarge). Reported on it are:

Curvature WFS

Set wfs.type="curvature". The only other mandatory parameters to define are wfs.l, wfs.nsubperring and wfs.lambda:

CWFS Subaperture map

Pyramid WFS

Briefly: The Pyramid WFS, in yao, is defined very similarly to the SH WFS. You have to define wfs.shnxsub (# of subapertures in pupil diameter), and the number of pixels in the pupil/phase plane per subaperture, wfs.npixpersub. For instance, if you want to simulate a pyramid WFS with 8x8 "subapertures" (pixels in pupil image), you can use: wfs.shnxsub=8, wfs.npixpersub=8, in which case sim.pupildiam=64. Currently, we can't have sim.pupildiam != wfs.shnxsub * wfs.npixpersub.

Other parameters used by the pyramid WFS only are listed in the structure table below. They are commented and fairly self-explanatory. Note that it is important (and required) to set a field stop when working with pyramid wfs in yao. The way it's implemented, the execution speed depends on the field stop diameter, so don't be surprised. Use a field stop large enough to avoid aliasing and growth of waffle like modes.

Zernike WFS

Set wfs.type="zernike". The only other mandatory parameter is wfs.nzer, which, as its name says, define the number of zernike in the WFS. Note that this include piston, so wfs.nzer=11 would include up to spherical (included). Important note: The Zernike are the Noll Zernike (see Noll 1976), which are not the standard Zernike, but the one used most commonly in AO.

User defined WFS

Set wfs.type="name_of_your_function".

Users can define their own WFS, to be integrated into the yao flow. The API are relatively straighforward:

func user_wfs_func(ipupil,phase,ns,init=)

The function should return a measurement vector (show/hide example).

Photometry and noise

Photometry and noise are only relevant (and implemented) for Shack-Hartmann and Curvature WFS.

All the photometry is controlled with the wfs.gsmag, wfs.skymag and wfs.zeropoint. wfs.zeropoint is defined in photons/second/full_aperture (incident on the telescope aperture, i.e. after crossing the atmosphere). Note that this depends on the telescope diameter (may be it was a bad choice to define it like that, but now it's done). It is simple enough to convert known photometric zero point to yao zeropoint. The number of photons available for WFSing is derived from this zeropoint and the above mentioned magnitude following regular equations (see also wfs.optthroughput).

yao doesn't have a concept of ADU. Or rather, I am assuming one electron/ADU across the board. So the signal in in electrons, and the noise is also in electrons.

Deformable Mirrors

The deformable mirrors are entirely and solely defined through their influence functions (shape and location). During the aoinit() phase, yao calls DM influence function definition functions, for each DM in the parameter file. Once computed, the influence functions are stored in the dm structure (*dm._def) and saved on disk for future runs (name-if#.fits files). Note that if you have a set of commands, you can get the DM surface by calling

surf = compute_dm_shape(nm,command_ptr);

where command_ptr is a pointer to a float vector containing the DM commands (e.g. command_ptr = &(float(indgen(dm._nact)))).

Stackarray (PZT,SAM,Piezo) DM

Set dm.type="stackarray". The only other mandatory parameters for stackarrays are dm.pitch and dm.nxact.


dm.nxact sets the number of actuators in the diameter of the DM (along the X or Y axis). This includes extrapolated/slaved actuators/guard rings.

dm.pitch sets the pitch (i.e. interactuator distance) in near field pixels.

Note that you have to dimension the stackarray DM according to your base system definition (sim.pupildiam). If you set sim.pupildiam=60 and dm.nxact=7 and you do not desire an extra ring of actuator, then you should set dm.pitch=10 to span the entire pupil diameter (7 actuators, 6 spaces between actuators...).

Influence functions

This was the first type of DM implemented in yao. With time, we have implemented several type of influence functions. By default, it uses a functional form which was first derived by J.-P. Gaffard of then CGE (now CILAS) to fit ZIGO measurements of CGE DM's real influence functions. It is possible to use newer (but not necessarily better) forms by setting dm.irexp:
dm.irexpfunctional form
0Old, adhoc functional form fitted on actual ZIGO data

irfact is set through dm.irfact.

Actuator Coupling

You can set the inter-actuator coupling with dm.coupling. dm.coupling=0.3 means 30% coupling between an actuator and its nearest neighbor. Typical value are 0.1-0.4. The almost universal consensus nowadays is to use coupling values of about 0.3.

Saving RAM & Disk with dm.elt

Stackarray mirrors have very local influence functions. Therefore, for large number of actuators, most of the surface of the DM is zero. It is thus a big waste of RAM and disk (these influence functions take a lot of space). Enabling dm.elt will save most of this space back: only a small subsection of the whole DM surface is kept, surrounding the actuator. The coordinates of this subimage is kept, and used later by comp_dm_shape() to compute the correct DM shape. The call to comp_dm_shape() is transparent. Depending on the value of dm.elt for the current DM, comp_dm_shape() will call the relevant C routine. In fact, there is no compromise in setting dm.elt, but its benefit really only show up for large stackarray DM (dm.nxact≥10).

Extrapolated/Slaved actuators

The idea in yao is to define all of your stackarray mirror in your parameter file, then run aoinit() and calibrate the iMat, then select which actuators are controlled and which are extrapolated based on the WFS sensitivity to it. Use dm.thresholdresp to select this threshold. This is a fractional number, so 1 would filter all actuators and 0 would retain them all. 0.5 would retain all actuator which measurement max response is larger than the max response of all actuators. Note that if you set dm.thresholdresp to a negative value, aoinit() will enter an interactive mode into which you can test various values of dm.thresholdresp (Show/Hide example).

Once extrapolated actuators are selected, yao stores the valid influence functions and the extrapolated ones in 2 different files (e.g. test-if1.fits and test-if1-ext.fits). Influence functions are not recomputed in further calls to aoinit() unless the clean keyword is set (in which case aoinit() restart from scratch and recompute everything).

Bimorph (curvature) DM

Set dm.type="bimorph". The other mandatory parameter is dm.nelperring.

dm.nelperring: This is similar to wfs.nsubperring (see here) and defines the number of electrodes per ring.

Now, most curvature mirror are a bit different from this simple design provided by defining dm.nelperring, in which there is no gap between rings. Many curvature DMs have for instance a large gap between the outer ring of electrodes (the one that create mostly boundary conditions) and the next inner electrode ring. To create more realistic DM electrode pattern, use angleoffset, rint, rout, supportRadius and supportOffset.

The general idea follows what has been exposed in the curvature WFS area above. With rint and rout, one can define the inner and outer radius of each ring of electrodes, expressed in fraction of the pupil radius. With angleoffset, one can define a global rotation per ring, and with the support* parameters, one can define the physical properties of the 3 DM support points.

That brings us to the next subject: how curvature DMs are modelled. Right now, what we do is solve the Poisson equation, imposing a given curvature over a given electrode and integrating twice to get the DM surface. Then, we impose that the three points defined as the support location go through the z=0 plane. It's that simple. It means that we do not currently support fancier support method like continuous support (e.g. by a rubber ring), 6 points support, etc..

Modal (Zernike, Disk Harmonic, Karhuenen-Loeve) DM

Set dm.type="zernike","diskharmonic" or "kl". Other mandatory keywords are dm.nzer (for zernike), dm.max_order for disk harmonic or dm.nkl (for KL), which specify the number of modes. Note that nzer (Zernike) & max_order (disk harmonic) will start and include piston, while nkl (KL) will not (there is formally no piston mode in the KL basis).

Segmented DM

Set dm.type="segmented". The other mandatory keywords is dm.nxseg.

dm.nxseg is the number of segment in the long diameter (X axis).

dm.fradius stands for "filter radius". Segments are created over a wider area than the nxseg defined above. Only segments which distance to the (0,0) pupil coordinates is ≤ fradius will be kept. default dm.pitch*dm.nxseg/2.

User defined DM

Set dm.type="name_of_your_function" (as a string).

The API for the user provided DM function are as follow:

func make_my_dm(nm,disp=)
where nm is the DM number. Use the disp keyword to display whatever you need to display when this routine is called. The goal of this routine is to create maps of the influence functions. These maps have to be stored (thus this function has to create and fill) in the array *dm._def (dm._def is a pointer to the data). The proper way to do that is:

func make_my_dm(nm)
  extern dm;   // dm is extern (global)
  dim = dm(nm)._n2-dm(nm)._n1+1; // <- important
  def = array(float,[3,dim,dim,n_actuator]);
  for (i=1;i<=n_actuator;i++) def(,,i) = some_form; // <- fill it to your taste
  dm(nm)._def = &def;  // <- store influence functions
  dm(nm)._nact = n_actuator;  // <- important, number of actuators

Explanation: In this user supplied function, you have to define the influence functions (*dm._def) and, importantly, the number of actuators (dm._nact). Prior to the call to this function, aoinit() has defined dm._n1 and dm._n2, which are the indices at which these influence function will go in the big sim._size x sim._size array used for FFT, etc. The influence functions are reduced in size to save RAM, as they can potentially take a lot of memory space. Normally, dm._n1 and dm._n2 are defined to include the pupil pixel diameter (sim.pupildiam) plus 4 pixels padding on each side if the DM altitude is zero (in the pupil), and are defined to span the entire sim._size array if the DM altitude is above ground (as there is no simple way to know where to stop exactly as this depends on the WFS and target locations). Note that dm is defined in extern as you are modifying its value.

Tip-Tilt Mirror

Set dm.type="tiptilt".

This creates a regular tip-tilt mirror. The units in the command vector are arcsecs.

Anisoplanatism "DM"

Set dm.type="aniso".

This has been implemented because of the need in MCAO to control plate scale modes (also called anisoplanatism modes or quadratic modes). You will need to set dm.alt to the altitude of an existing DM of a regular type (e.g. stackarray) that will be responsible for the creation of the anisoplanatism modes (see mcao2-bench.par for an example of aniso DM).

Hysteresis and misregistration/offset


Hysteresis is implemented at the top level so that it works for any kind of DM (during the AO loop only). The parameter dm.hysteresis is fairly self explanatory: dm.hysteresis=0.1 means 10% hysteresis.


dm.misreg allow to induce DM misregistration (on the fly). dm.misreg is expressed in pixels, so you will have to do a small conversion to get it in other units, but this way, it works for any kind of DM. Note that dm.misreg can be fractional (fraction of a pixel). Because the size of the array containing the DM influence function is limited to save memory space, the misregistration using this parameter can only be of a few pixels (up to 4 I believe for a ground-conjugated DM). If you wish to induce more registration, use the parameter dm.pupoffset (see below).


Yao v4.5.0 introduced a new kind of influence function offset parameter: dm.pupoffset. dm.pupoffset is entered in meters (a 2 elements vector which contains the x and the y offsets). Those are transversal offset of the influence functions in the near field. This functionality is implemented in the compute_dm_shape() function, so it is general and works for all types of DM. dm.pupoffset is intended for large offsets, and the influence functions are shifted by an integer number of pixels (the requested offset in meter, rounded to the nearest pixel, see sim.pupildiam and tel.diam). If you want finer control over the DM position, use dm.misreg (see above). dm.pupoffset is similar to wfs.pupoffset (see WFS section), but for the DMs.

Building a Reconstructor

The default reconstructor which was initially the only option in YAO is the truncated singular value decomposition (SVD). This is a least-square solution where invisible or poorly sensed modes are totally discarded. A faster and better performing reconstructor (especially for large AO systems) is the minimum mean-square error (MMSE) reconstructor. This MMSE reconstructor can be either built using full matrices or sparse matrices techniques, which offers a significant speed boost.

Method selection

There are three options for wavefront reconstruction, and these are set using the mat.method parameter.

Truncated singular value decomposition (default)
mat.method = "svd"
mat.condition = &([desired conditioning])
Minimum Mean Square Error Reconstructor

This method computes the reconstructor as

R = (HT H + a C)-1 HT

where H is the interaction matrix (poke matrix), a is a regularization parameter and C is a regularization matrix. The regularization matrix (see dm.regtype or dm.regmatrix) can be user specified or computed by YAO. If it is computed, it will use an identity matrix by default, or a matrix created by convolving a laplacian operator by itself. This matrix has similar statistics to the inverse covariance matrix for Kolmogorov turbulence and penalizes local waffle in the DM.

The relevant parameters for the MMSE inversion process are

    mat.method = "mmse";
    dm.regparam = 0.01; // higher for larger systems
    dm.regtype = "laplacian"; // "laplacian" or "identity"
    dm.regmatrix = &([]); // you can specify your own regularization matrix
Minimum Mean Square Error Sparse Reconstructor

A sparse matrix implementation of the above is also available. This is ideal for very large systems with thousands of actuators, as it reduces memory requirements and increases computational speed. Do not use the sparse method when the interaction matrices are not actually sparse (e.g., when using a Zernike DM).

  mat.method = "mmse-sparse";
  mat.sparse_MR = 1e4; // a number larger than the total number of rows
  mat.sparse_MR = 1e5; // a number larger than the number of elements
  dm.regparam = 0.01; // higher for larger systems
  dm.regtype = "laplacian"; // "laplacian" or "identity"
  dm.regmatrix = &([]); // you can specify your own regularization matrix

If you set the values of mat.sparse_MN and mat.sparse_MR too low, the system could crash. Too high and you could run out of memory.

Code Your Own Reconstructor

It is possible to code your own reconstructor and make YAO use it. You need to overwrite the reconstructor matrix cMat as the matrix that you compute. You can also save the reconstructor so that it is not overwritten whenever aoinit is run.

    cMat = myNewReconstructor;

Your reconstructor will, however, be overwritten if you do aoinit,clean=1 or aoinit,svd=1

Loop parameters

Most AO systems run in closed-loop. Typical parameters are
  loop.method     = "closed-loop";
  loop.framedelay = 1;
  loop.gain       = 0.6;
  loop.leak       = 0.001; 

The gain is typically a value around 0.5, and the leak is a small value needed to prevent the build-up of poorly sensed modes on the DM. The framedelay is the compute time rounded to the nearest number of frames, typically 1. There is already a one cycle delay inherent in the zero-order hold of the WFS and the DM stare, and this does not need to be added.

For tomographic reconstructors, you can use pseudo open-loop control: loop.method = "pseudo open-loop". If you have an AO system where the WFSs see the uncorrected wavefront, you can run the simulation in open-loop:

    loop.method     = "open-loop";
    loop.framedelay = 1;
    loop.gain       = 1;
    loop.leak       = 1.; 

In open-loop, it is essential to set the loop gain to 1 and the leak to 1.

It is also possible to create high-order control laws and to control each DM with a separate controller. This is especially useful to separate the tip-tilt and DM control. For example, to implement the control law on dm(1):

$$ \frac{0.5 + 0.12 z^{-1} -0.05 z^{-2}}{1 - 0.64 z^{-1} -0.36 z^{-2}} $$
dm(1).ctrlnum = &([0.5,0.12,-0.05]);
dm(1).ctrlden = &([1,-0.64,-0.36]);

Running yao using parallel processes

You can take advantage of multi-core machines to speed up the execution time. Gains range from 20% to 350% (i.e. 1.2x to 3.5x). Here is a copy of the blog entry I posted when announcing yao parallel:

More work on yao parallel:

  • I have merged some of the -useful- parallel features implemented in our Gemini local yao_mcao to the yao trunk. Namely,
    • The DM/WFS parallelization (meaning, taking advantage of the fact that real systems have at least one frame delay (remember that even if the command is applied immediately after the measurements are received, there is still one frame delay due to the zero order hold). The one frame delay means the WFSing at iteration N is not needed to compute the command at iteration N, thus can be computed in parallel.
    • The PSF calculations, which is a no brainer.
  • I have implemented the parallelization of individual SHWFS. Meaning, an individual SHWFS can now run in parallel on an arbitrary number of "threads" (in fact, forks). This works only for Shack-Hartmann WFS. The good thing is that this is implemented fairly deep into the code, so any call to sh_wfs() will benefit from this, which means the iMat acquisition is accelerated too.
Depending on your system, i.e. WFS order, number of WFSs/DMs, whether you use dm.elt or not, sparse reconstructor or not, and more importantly, on the machine you run yao on, using the parallel features can make you gain anywhere between 20% and 200% in execution speed. If you want to use it, you will most probably have to experiment a bit to see what option should be set. Guidelines are given in the manual. The bottom line is that you should not use more forks (processes) than you have CPUs, otherwise swapping processes in/out of L2 cache introduce an overhead which cancels the paralleization advantage. With my core2duo laptop, I can typically gain 30 to 50%. With an 8 cores Xeon-based machine we have at work, I gain typically close to a factor of 2 (can be more or less). This does not seem particularly impressive, but one has to consider that a yao loop is made of many different individual operations: WFS, cMat product to retrieve the error, DM shape computation from the DM commands and the influence functions, PSF calculation, etc... So a gain in one of this items is only going to benefit the whole thing marginally. As far as the individual WFS parallelization is concerned, here are some hard numbers:
8 Core Xeon (2 boards) @ 3.33GHz, running RHEL 5, 32 bits. yorick/yao built with -O2.
16x16 SHWFS system, physical model, no noise. All times (avg & median) in milliseconds for one frame.

sh16x16_svipc:        1 threads, sh_wfs() avg=2.86, median=2.90
sh16x16_svipc:        2 threads, sh_wfs() avg=1.87, median=1.89
sh16x16_svipc:        4 threads, sh_wfs() avg=1.37, median=1.37
sh16x16_svipc:        8 threads, sh_wfs() avg=1.14, median=1.15
sh16x16_svipc:       16 threads, sh_wfs() avg=2.33, median=2.30
You can see that the speed increases up to 8 forks (max gain = 2.5x). Adding more forks makes the thing slower. The gain does not reach larger value just because (1) there are overheads associated to the use of shared memory/semaphores and (2) there is also a global overhead in the sh_wfs() function. What's parallelize is the core of sh_wfs() (FFTs, noise application and slopes calculations from spots), but any overhead within the sh_wfs() function itself, before calling the FFT/slopes calculation routine is not going to gain from the parallelization.

The gain brought by the other, more general parallel features (DM/WFS & PSFs) is nicely illustrated by the following examples. This was run on the same machine (8 cores Xeon), with curvature.par (36 subapertures curvature system ran on a 120 pixels pupil). Look at the following results (I have omitted the timing results for items that are not of concerned here and were negligible):

sim.svipc=0 (no parallelization)
time(1-2) =  9.89 ms  (WF sensing)
time(3-4) =  0.01 ms  (cMat multiplication)
time(4-5) =  1.04 ms  (DM shape computation)
time(5-6) =  5.84 ms  (Target PSFs estimation)
59.18 iterations/second in average

sim.svipc=1 (WFS/DM parallelized)
time(1-2) =  4.11 ms  (WF sensing)
time(3-4) =  0.01 ms  (cMat multiplication)
time(4-5) =  0.94 ms  (DM shape computation)
time(5-6) =  5.78 ms  (Target PSFs estimation)
91.43 iterations/second in average

sim.svipc=2 (PSFs parallelized)
time(1-2) = 10.69 ms  (WF sensing)
time(3-4) =  0.01 ms  (cMat multiplication)
time(4-5) =  1.05 ms  (DM shape computation)
time(5-6) =  1.23 ms  (Target PSFs estimation)
76.41 iterations/second in average

That's a typical example. The "no parallel" case is dominated by WFSing (10ms) and PSF (6ms) calculations. When parallelizing the WFS (sim.svipc=1), the WFSing time goes down to 4ms. Indeed, now 6ms of it can be done in parallel with other tasks: here, the PSF calculation that takes 6ms. With sim.svipc=2 (only the PSF are parallelized), the PSF calculation time goes to basically to 0 while the WFSing stays at 10ms. It all makes sense. I haven't put the sim.svipc=3 results as they are a little more difficult to understand and would just introduce confusion.

The individual WFS parallel feature is controlled by the parameter wfs.svipc. The "global" DMs/WFSs abd PSF parallelization is controlled by sim.svipc (bit 0 controls DM/WFS and bit 1 controls PSFs, thus sim.svipc=3 means both are turned on).

Of course, if you want to use yao parallel facilities, you will need to install the yorick-svipc plugin, through most of the normal channels. It runs on Linux (extensively tested) and OsX (somewhat tested).

As it says, you need to install yorick-svipc. It is available right now (July 2010) through the maumae.net package source web page and on archlinux AUR. It will eventually make it in the Debian & Ubuntu repos.

Parallel implementation gotcha & bugs

It is difficult to code for parallel deployment. Especially for something like yao which offers a lot of various configuration option to the users. I have tested the parallel implementation quite a lot, but I am sure I have left some corners unexplorered/untested. So you might hit some unexpected features. Let me know if you do notice something weird. In the meantime, here is a list of things to watch out for when using the parallel (svipc) thingy (the first three items should not happen. I haven't seen it for a long time. But I originally saw it hence I give you these tricks in the unlikely event it happens) :

Other Features

More to come (a note to myself).

Yao Structures

Comment on notation: &float = pointer to a float array/vector (idem for &long, &string, etc...). Define it as:
atm(1).layerfrac = &([0.5,0.2,0.3]);
I know, it seems weird, but there is a good reason to have done it like this.
sim structure
name string N/A none no A name for this simulation run
pupildiam long pixels none yes Pupil diameter
debug long N/A 0 no Debug level
verbose long N/A 0 no Verbose level
svipc long N/A 0 no bitwise set features for parallelization:
0 = no parallelization
bits effect
0(1) WFS/DM global split (2 process)
1(2) PSFs calculation parallelization
e.g. sim.svipc = 1 -> split WFS/DM
e.g. sim.svipc = 2 -> parallelize PSFs
e.g. sim.svipc = 3 -> WFS/DM & PSFs
shmkey long N/A no shared memory key. Change in each session if running multiple yao session simultaneously.
semkey long N/A no semaphore svipc key. Change in each session if running multiple yao session simultaneously.
atm structure
dr0at05mic float Unitless none yes Dr0 at 0.5 microns, at zenith
screen &string N/A none yes Phase screen file names
layerfrac &float Unitless none yes Layer fraction. Sum to one is insured in aoinit
layerspeed &float meter/sec none yes Layer speed
layeralt &float meter none yes Layer altitude, at Zenith
winddir &long Unitless 0 yes Wind dir (not operational, use 0 for now)
wfs structure
type string N/A none yes Valid types are "curvature", "hartmann", "pyramid", "zernike" or "user_function" where user_function is the name of a function defined by the user (see doc)
lambda float micron none yes WFS wavelength in microns
subsystem long N/A 1 no Subsystem this WFS belongs to
gsmag float Unitless 0 no WFS guide star magnitude. For LGSs, see below. Beware to set to something bright when doing your iMat
skymag float Unitless 0 (no sky) no WFS sky magnitude/arcsec^2. [0 means no sky]
noise long N/A 0 no Enable noise (photon noise & read out noise)
ron float e- 0 no Read out noise
darkcurrent float e-/s/pixel 0 no Dark current
gspos(2) float arcsec [0,0] no WFS guide star position [x,y]
gsalt float meter 0 no WFS guide star altitude, at zenith. 0 for infinity.
gsdepth float meter 0 no WFS guide star depth in meter, at zenith (e.g. Na layer thickness)
laserpower float Watt 0 See commentWFS laser power projected on sky (Na laser only). Required when using lasers. Exclusive with gsmag i.e. define one OR the other
filtertilt long N/A 0 no Filter TT on this sensor (0=no)
correctUpTT long N/A 0 no Correct uplink tip-tilt (0=no)
uplinkgain float Unitless 0 no Uplink TT loop gain
dispzoom float N/A 1.0 no Zoom factor for the display, useful in multi-WFS configuration (typically around 1)
optthroughput float Unitless 1.0 no Optical throughput to WFS
disjointpup long N/A 0 no If set, the WFS#n will be filtered by an array disjointpup(,,n) that has to be defined by the user. see user_pupil(). Allow for GMT-type topology
Pyramid WFS only keywords
pyr_mod_ampl float arcsec 0.0 yes pyramid wfs modulation amplitude radius [arcsec]
pyr_mod_npts long N/A 1 yes total number of point along modulation circle [unitless]
pyr_padding long N/A 0 no Pad the pupil image to reduce spatial aliasing [unitless]. A pad of 1 means adding wfs.npixpersub pixels on each side of the pupil image. Typical 0 to 4.
pyr_mod_loc string N/A "after" yes Location of modulation (before/after the field stop). Valid values are "before" or "after"
Curvature WFS only keywords
l float meter none yes Extra focal distance in a F/60 beam
nsubperring &long Unitless none yes # of subapertures per ring. See doc
angleoffset &float degree 0 no CCW Offset angle for first subaperture of ring. See doc
rint &float Unitless See commentno Inner radius for each ring, in fraction of pupil radius
rout &float Unitless See commentno Outer radius for each ring, in fraction of pupil radius
fieldstopdiam float arcsec 1.0 no Diameter of field stop. Used only to compute sky contribution (with skymag)
Shack-Hartmann WFS only keywords
shmethod long N/A none yes 1= Geometric, simple gradient average over subaperture. 2=Diffractive, full propagation
shnxsub long Unitless none yes # of subapertures in telescope diameter
pixsize float arcsec none yes Focal plane: Subaperture CCD pixel size in arcsec
npixels int Unitless none yes Focal plane: Final # of pixels per subaperture
npixpersub long Unitless none no Pupil plane: # of pixel in a subaperture (to force npixpersub and bypass constraint that pupildiam should be a multiple of this number e.g. to investigate lenslet larger than pupildiam
pupoffset(2) float meter [0,0] no Pupil plane: Offset of the whole WFS subapertures w.r.t telescope aperture. Allow misregistration w.r.t telescope pupil and other funky configurations
shthreshold float e- 0 no Threshold for the computation of the subaperture signal from CCD spots >= 0
biasrmserror float e- 0 no rms error on WFS CCD bias in electron
flatrmserror float Unitless 0 no rms error on WFS CCD flat, referenced to 1 (i.e. 0.1 mean 10% rms error). Typical value can be 0.01
fsname string N/A none no Fits file with subaperture amplitude mask. It should have dimension 2^sdimpow2 square. Can be float or long.
fstop string N/A "square" no Valid fiels stop type are "none", "square" or "round"
fssize float arcsec sub. size no Side (square) or diameter (round) of field stop
fsoffset(2) float arcsec [0,0] no Field stop offsets [x,y]
kernel float arcsec See commentno FWHM in arcsec of WFS spot gaussian kernel. Default is computed as a function of D/r0 and only used during iMat calibration
nintegcycles int Unitless 1 no # of iterations over which to integrate before delivering slopes
fracIllum float Unitless 0.5 no Focal plane: Fraction of subaperture illuminated for the subaperture to be valid
LLTxy(2) float meter [0,0] no Coordinates [x,y] of the laser projector, if any
centGainOpt long N/A 0 no Centroid gain optimization flag. Only for LGS (correctupTT and filtertilt must also be set for this to work)
rayleighflag int N/A 0 no Take rayleigh into account?
svipc long N/A 0 no Number of parallel processes (forks) to use for this WFS (0 or 1: don't parallelize, N: use main + (N-1) forks)
Zernike WFS only keywords
nzer int Unitless none yes Number of Zernike to be sensed. Starts at piston included.
dm structure
type string N/A none yes Valid types are "bimorph", "stackarray" "tiptilt", "zernike", "diskharmonic", "kl", "segmented", "aniso" or "user_function", where user_function is the name of a function provided by the user
pitch long pixel none yes Actuator pitch (pixel). stackarray/segmented only
subsystem long N/A 1 no Subsystem this DM belongs to
iffile string N/A "" no Influence function file name. Leave it alone.
alt float meter 0 no Conjugation altitude, at zenith
hyst float Unitless 0 no DM actuator hysteresis (0. to 1.)
push4imat float Volt 20 no Voltage to apply for iMat calibration. Note: the default is not OK for many configs. Change at will
thresholdresp float Unitless 0.3 no Normalized response threshold (in WFS signal) below which an actuator will not be kept as valid
unitpervolt float mic/Volt 0.01 no Influence function sensitivity in unit/volt. Stackarray: micron/Volt, Tip-tilt: arcsec/Volt.
maxvolt float Volt none no Saturation voltage (- and +) in Volt. None if not set
gain float Unitless 1.0 no Loop gain for this DM (total = this times loop.gain)
misreg(2) float pixel [0,0] no DM misregistration [x,y]
xflip long N/A 0 no Flip DM left/right (0=no)
yflip long N/A 0 no Flip DM up/down (0=no)
pupoffset(2) float meter [0,0] no Global offset of whole actuator pattern w.r.t pupil
disjointpup long M/A 0 no If set, dm(n) will be filtered by an array disjointpup(,,n) that has to be defined by the user. See user_pupil(). Allow for GMT-type topology.
Stackarray-only (SAM, PZT) keywords
nxact long Unitless none yes Number of actuator in pupil diameter
elt long N/A 0 no ELT mode: allow to save huge amount of RAM and time for the computation of the DM shape. No drawback
coupling float Unitless 0.2 no Influence function coupling coefficient
ecmatfile string N/A none no Valid to extrapolated projection matrix (extrap_com)
noextrap long N/A 0 no Set to disable use of extrapolated actuators
pitchMargin float Unitless 1.44 no Margin to include more corner actuators when creating inf.functions optional [1.44]
irexp long N/A 0 no Use original functional form (irexp=0) or exp(-(d/irfact)^1.5) model (irexp=1) or sinc*gaussian (irexp=2)
irfact float Unitless 1.0 no use when irexp=1 (see above)
regparam float Unitless 0.0 no Parameter that determines the strength of the regularization. Set to zero for systems like tip-tilt mirrors and to a small positive value (e.g., 1e-4) for stackarray DMs. Experiment by trial an error to get the value that produces the highest Strehl. Experience shows that any value within a factor of four of optimal produces the same results.
regtype string N/A "identity" no Regularization approach: "identity" (default) or "laplacian"
regmatrix &float Unitless none no Pointer to the regularization matrix to use instead of computing one of the two options. This is usually not necessary, but is left as an option for advanced users. A good choice would be the inverse covariance matrix of the turbulence, for example.
Bimorph-only keywords
nelperring &long Unitless none yes Number of electrodes per ring, e.g &([6,12,18])
angleoffset &float degree 0 no Offset angle for first electrode in ring
rint &float Unitless See commentno Inner radius for each ring, see doc
rout &float Unitless See commentno Outer radius for each ring, see doc
supportRadius float Unitless 2.2 no Radius of DM 3 support points, normalized in pupil radius
supportOffset float dgree 90 no Angle offset of first support point
Zernike DM-only keywords
nzer long Unitless none yes Number of Zernike modes, including piston
Disk Marmonic-only keywords
max_order long Unitless none yes Number of Disk Harmonic modes, including piston
Karhuenen-Loeve DM-only keywords
nkl long Unitless none yes Number of Karhuenen-Loeve modes, including piston???
Segmented DM-only keywords
nxseg long Unitless none yes Number of segments in long axis (X)
fradius float pixel See commentno Segments are created over a wider area than the nxseg defined above. Only segments which distance to the (0,0) pupil coordinates is <= fradius will be kept. default dm.pitch*dm.nxseg/2.
mat structure
method string N/A "svd" no Reconstructor creation method ("svd", "mmse" or "mmse-sparse")
condition &float Unitless none yes Condition numbers for SVD, per subsystem.
file string N/A none ??? iMat and cMat filename. Leave it alone.
sparse_MR long Unitless 10000 no Sparse only: maximum number of rows (actuators) in the imat
sparse_MN long Unitless 200000 no Sparse only: maximum number of elements in the imat
tel structure
diam float meter none yes Telescope diameter
cobs float Unitless 0. no Central obstruction / telescope diameter ratio
Tip vibration parameters
tipvib_white_rms float arcsec 0. no rms of vibration white noise
tipvib_1overf_rms float arcsec 0. no rms of vibration 1/f noise (from 1 Hz to cutoff)
tipvib_peaks &float Hz 0 no positions of vibration peak in PSD
tipvib_peaks_rms &float arcsec 0 no rms of each vibration peaks (defined above)
tipvib_peaks_width&float Hz 1freq binno width of each vibration peaks (default 1 freq bin)
tilt vibration parameters
tiltvib_white_rms float arcsec 0. no rms of vibration white noise
tiltvib_1overf_rmsfloat arcsec 0. no rms of vibration 1/f noise (from 1 Hz to cutoff)
tiltvib_peaks &float Hz 0 no positions of vibration peak in PSD
tiltvib_peaks_rms &float arcsec 0 no rms of each vibration peaks (defined above)
tiltvib_peaks_width&float Hz 1freq binno width of each vibration peaks (default 1 freq bin)
target structure
lambda &float micron none yes Image wavelengths in micron
xposition &float arcsec none yes "Targets" X positions in the field of view
yposition &float arcsec none yes "Targets" Y positions in the field of view
dispzoom &float Unitless 1. no Display zoom, useful for multi-targets. Typically around 1
gs structure
zeropoint float See commentnone yes Photometric zero point (#photons@pupil/s/full_aper, mag0 star).
zenithangle float degree 0 no Zenith angle. The zenith angle is used to compute: r0 off-zenith, atmopheric turbulence layer altitude, LGS altitude and thickness of Na Layer, LGS brighness note that dm altitude is unchanged
lgsreturnperwatt float phot/cm2/s 22. no Sodium LGS return in photons/Watt/cm2/s at entrance pupil, at zenith. Modified by gs.zenithangle. Basically, you have to fold in this the sodium density and your model of return
loop structure
niter long Unitless none yes Total number of iterations
ittime float second none yes Iteration time (sampling time)
gain float Unitless 0 See commentLoop gain. Optional, but important!
leak float Unitless 0 no Leak term (0 means no leak)
gainho &float Unitless 0 no Higher order gains (starting at 2nd order, up to 10th)
leakho &float Unitless 0 no Higher order leaks (starting at 2nd order, up to 10th)
framedelay long frames 0 no Loop delay (# of frames). Regular CCD 1 frame integration -> framedelay=1 + readout & Calculation -> framedelay=2
startskip long iteration 10 no Number of iteration to skip before collecting statistics
skipevery long iteration 0 no In phase screen, skip by "skipby" steps every "skipevery" iterations (0=none). See doc
skipby long iteration 10000 no See above. This is to get better statistical coverage
stats_every long iteration 4 no Compute stats every so many iteration
jumps2swapscreen long Unitless 0 no Number of jumps (i.e. niter/skipevery) after which screens will be swapped (rotation, 2->1, 3->2... 1->last). Default is no jump.
modalgainfile string N/A "" no Name of file with mode gains
opt structure
phasemaps string N/A none no Filename of phasemap. Z scale should be in microns
alt float meter 0 no Equivalent altitude in m.
misreg(2) float pixels [0,0] no Misregistration [x,y] (similar to DM, see above)

Scripting and Hacking Yao


Thanks to yorick and the structure of yao, scripting is relatively easy. The principle is as follow:

Here is a simple example you can find in the examples directory (yao_loop_example.i):


// check if generic phase screens exist, if not, create them:
if (!open(Y_USER+"data/screen1.fits","r",1)) {

// the "wait" is needed here:

// read out parfile
atm.dr0at05mic = 35; // be more gentle

// Define vector on which we want to loop and final strehl array.
// We want to estimate performance for 3 values of the guide star
// magnitude and 4 values of the loop gain (for instance)
gsmagv = [6,9,12];           // guide star mag vector
gainv  = [0.01,0.1,0.5,1.0]; // gain vector
// strehl array to store results:
strehlarray = array(0.,[2,numberof(gsmagv),numberof(gainv)]);

// loop on gsmag and gain
for (ii=1;ii<=numberof(gsmagv);ii++) {
  for (jj=1;jj<=numberof(gainv);jj++) {
    // it's safer, but not always necessary, to call again
    // aoinit (here for gsmag). some parameters do not need it.
    // Setup loop:
    // go: do all loop.niter
    go, all=1;
    // after_loop() is now called automatically at last iter of go()
    strehlarray(ii,jj) = strehllp(0); // fill in result array
    // and display results as we go:
    for (ll=1;ll<=ii;ll++) {
      ylims=limits()(3:4); ymax=ylims(2); yspace=(ylims(2)-ylims(1))/15.;
      plt,swrite(format="gsmag=%d",gsmagv(ll)),0.011,ymax-yspace*(ll-1), \
    xytitles,"Loop Gain",swrite(format="Strehl @ %.2fmicrons",(*target.lambda)(0));

You can do all kind of scripting like this to find for instance the optimal working point for a system in given conditions (above, what is the best gain for a given GS magnitude). We could also have added a loop on wfs.ittime in the example above. One has however to be careful: Probing multi-dimensional spaces can quickly be overwhelming in execution time. Generally, I consider that 10000 iterations are needed to give a statistically significant answer. So for instance, in the example above. that would mean 4x3x10000=120000 iterations total. At 90 iterations/seconds, this means 1300 seconds (about 20mn). Adding 4 points for wfs.ittime would lead to 480000 iterations, or 1 hour and 20mn. This can quickly become prohibitive.

Hacking Yao

At one point or another, if you are serious about simulating your system, it is likely that you will need this or that feature, or modification to the existing yao code.

Go ahead and dive. Feel free to modify yao. It's open source after all. The goal of this too-short section is to give you some head start for doing just that: hacking yao.

So, various points, in no particular order:


There is a lot of things missing in this documentation. I will try to complement it as time allows. I hope however that this will help you getting the best out of yao and avoiding frustration. Yao is a fairly complete AO simulation tool. It is flexible and fast. Flexibility means I can not exert too much control over the user input parameters. If I were to do that, I would necessarily impose a straitjacket over your creativity, and we want to avoid this at all cost. So, yes, there is a steep learning curve, but at the end, you should be able to wrestle yao into submission! I have only one more advice to give: check your results. And recheck them. And make sure it makes physical sense.