.. _scao_basic_tutorial:

SCAO Basic Tutorial: Complete Walkthrough
=========================================

This comprehensive tutorial guides you through creating and running a basic Single Conjugate Adaptive Optics (SCAO) simulation using SPECULA.
A more advanced SCAO simulation tutorial is available in the :ref:`scao_tutorial`.

**What you'll learn:**

* Setting up a minimal SCAO system configuration with a Pyramid WFS
* Running calibration and closed-loop phases
* Using Zernike modes and a single atmospheric layer
* Setting the push-pull amplitude directly in the YAML file

**Prerequisites:**

* SPECULA installed and working (see :doc:`../installation`)
* Basic understanding of adaptive optics concepts
* Python and YAML familiarity

.. note::

   **About atmospheric and source parameters:**
   All atmospheric parameters (seeing, layer heights) and source heights are defined at zenith.
   See :ref:`simulation_parameters` for details on units, conventions, and the zenith angle treatment.


Tutorial Overview
-----------------

We'll simulate a simple SCAO system with:

* Circular pupil (VLT-like, 8m class)
* 1 atmospheric layer
* Pyramid wavefront sensor
* Zernike modal DM
* 1 kHz control loop with integrator controller
* R-band natural guide star (magnitude 8)

Part 1: System Configuration
----------------------------

In this part, we will define the system configuration in a YAML file.
SPECULA uses a modular, object-oriented architecture where each component of the simulation is defined as an object in a YAML configuration file.
In the YAML file all the parameters of the simulation are defined, including the objects and their connections.

Understanding the YAML Simulation File
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The core of a SPECULA simulation is defined in a YAML file, which describes the system as a set of interconnected objects.
Each **section** in the YAML file corresponds to a specific object (or "block") in the simulation, such as the atmosphere, the wavefront sensor, the deformable mirror, or the control law.

**How it works:**

- **Object Instantiation:**  
  Each section's ``class`` parameter specifies which Python class will be instantiated for that object.
  All other parameters in the section are passed as arguments to the class constructor.
  This is handled automatically by the simulation engine (``simul.py``).

- **Inputs and Outputs:**  
  The ``inputs`` and ``outputs`` fields define how objects are connected.  
  
  - The value of each input is a reference to the output of another object, using the syntax ``object_name.output_name``.
  - This mechanism defines the **data flow** and the **dependency graph** of the simulation.

- **Simulation Graph:**  
  The entire simulation is built as a directed graph, where nodes are objects and edges are the connections defined by the ``inputs``.
  The simulation engine parses the YAML, instantiates all objects, and connects them according to this graph.
  See :doc:`../simul_diagrams` for visual examples of simulation graphs automatically generated by SPECULA.

- **Special Notes:**  
  
  - In the propagation block (``prop``), the input ``common_layer_list`` often includes ``dm.out_layer:-1``. This is a SPECULA convention to handle feedback in the simulation loop. The ``-1`` index is used to resolve an ambiguity in closed-loop simulations: not all the elements of the loop can perform their operations at the same time, at least one must happen at the following time step. In an Adaptive Optics context the DM output is computed during the current time step, but it is applied in the next time step.

**Example:**

.. code-block:: yaml

   prop:
     class: 'AtmoPropagation'
     simul_params_ref: 'main'
     source_dict_ref: ['on_axis_source']
     inputs:
       atmo_layer_list: ['atmo.layer_list']
       common_layer_list: ['pupilstop', 'dm.out_layer:-1']
     outputs: ['out_on_axis_source_ef']

In this example:

- The ``prop`` object is an instance of the ``AtmoPropagation`` class.
- It receives as input the atmospheric layers from ``atmo.layer_list``, the pupil geometry from ``pupilstop``, and the most recent DM surface from ``dm.out_layer:-1``.
- It produces an output ``out_on_axis_source_ef``, which can be used as input by other objects (e.g., the Pyramid WFS).

This modular and explicit approach makes it easy to customize, extend, and debug your simulation setup.

Create the main YAML configuration file
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Create a YAML configuration file, for example ``params_scao_pyr_basic.yml``:

.. code-block:: yaml

   # main section with simulation parameters used by most of the components
   # Note: zenith angle is not specified, so it is assumed to be 0 (on-axis) 
   main:
     class:             'SimulParams'
     root_dir:          './calib/'             # Root directory for calibration manager  
     pixel_pupil:       160                    # Linear dimension of pupil phase array
     pixel_pitch:       0.05                   # [m] Pitch of the pupil phase array
     total_time:        0.050                  # [s] Total simulation running time
     time_step:         0.001                  # [s] Simulation time step

   # Atmospheric seeing (approx. inverse of Fried parameter r0). Here is a static value,
   # but it can be a function of time. The full list of parameters can be found in the
   # init method of the WaveGenerator class.
   seeing:
     class:             'WaveGenerator'
     constant:          0.8                   # ["] seeing value (500nm and at zenith)

   # Wind speed and direction, also static in this example.
   # These can be functions of time or vary per layer in more complex setups.
   wind_speed:
     class:             'WaveGenerator'
     constant:          [20.]                 # [m/s] Wind speed value

   wind_direction:
     class:             'WaveGenerator'
     constant:          [0.]                  # [degrees] Wind direction value

   # Source definition. This defines the direction of the propagation the electromagnetic field
   # the flux intensity and the wavelength. In this case, we have a single on-axis source.
   # The source is defined in polar coordinates (arcsec, degrees). The full list of parameters
   # can be found in the init method of the Source class.
   on_axis_source:
     class:             'Source'
     polar_coordinates:  [0.0, 0.0]           # [arcsec, degrees] source polar coordinates
     magnitude:         8                     # source magnitude
     wavelengthInNm:    750                   # [nm] wavelength

   # Pupil stop definition. This is the pupil geometry used in our simulation.
   # It can be a circular aperture or a more complex shape.
   # It can be stored in a file or defined directly in the YAML.
   # When no parameters are specified, a circular pupil is assumed, with the size
   # defined by the pixel_pupil parameter in the main section. The full list of parameters
   # can be found in the init method of the Pupilstop class.
   pupilstop:
     class:             'Pupilstop'
     simul_params_ref:  'main'

   # Atmospheric layers generation and temporal evolution.
   # Here we define a single atmospheric layer with a constant Cn2 value.
   # Outer scale can be a scalar or a vector (for multiple layers).
   # The layer heights are defined in meters, and the Cn2 values must sum to 1.
   # The fov parameter defines the field-of-view in arcseconds.
   # The inputs are the seeing, wind speed, and wind direction defined above.
   atmo:
     class:                'AtmoEvolution'
     simul_params_ref:     'main'
     L0:                   40                   # [m] Outer scale
     heights:              [0.]                 # [m] layer heights at 0 zenith angle
     Cn2:                  [1.00]               # Cn2 weights (total must be eq 1)
     fov:                  0.0                  # [arcsec] Field of view, 0.0 means only on-axis,
                                                # no off-axis sources
     inputs:
       seeing: 'seeing.output'
       wind_speed: 'wind_speed.output'
       wind_direction: 'wind_direction.output'
     outputs: ['layer_list']

   # The propagation block simulates the propagation of the electromagnetic field
   # through the atmosphere, the pupil stop and the DM. It takes the source and the atmospheric layers
   # as inputs and outputs the electric field at the pupil plane in all the directions corresponding
   # to the source polar coordinates.
   # The output is a list of electric fields, one for each source direction.
   prop:
     class:                'AtmoPropagation'
     simul_params_ref:     'main'
     source_dict_ref:      ['on_axis_source']
     inputs:
       atmo_layer_list: ['atmo.layer_list']
       common_layer_list: ['pupilstop', 'dm.out_layer:-1']
     outputs: ['out_on_axis_source_ef']

   # The Pyramid WFS block simulates the Pyramid wavefront sensor.
   # It takes the electric field from the propagation block and computes the intensity on the detector.
   # The full list of parameters can be found in the init method of the ModulatedPyramid class.
   pyramid:
     class:             'ModulatedPyramid'
     simul_params_ref:  'main'
     pup_diam:          30.                     # Pupil diameter in subaps.
     pup_dist:          36.                     # Separation between pupil centers in subaps.
     fov:               2.0                     # Requested field-of-view [arcsec]
     mod_amp:           3.0                     # Modulation radius (in lambda/D units)
     output_resolution: 80                      # Output sampling [usually corresponding to CCD pixels]
     wavelengthInNm:    750                     # [nm] Pyramid wavelength
     inputs:
       in_ef: 'prop.out_on_axis_source_ef'

   # The detector simulates the CCD sensor where the Pyramid WFS intensity is recorded.
   # Its integration time can be a multiple of the simulation time step.
   detector:
     class:             'CCD'
     simul_params_ref:  'main'
     size:              [80,80]                 # Detector size in pixels
     dt:                0.001                   # [s] Detector integration time
     bandw:             300                     # [nm] Sensor bandwidth
     photon_noise:      True                    # activate photon noise
     readout_noise:     True                    # activate readout noise
     readout_level:     1.0                     # readout noise in [e-/pix/frame]
     quantum_eff:       0.32                    # quantum efficiency * total transmission
     inputs:
       in_i: 'pyramid.out_i'

   # The slope computer calculates the wavefront slopes from the detector frame.
   # It requires a pupil data object with the list of the valid sub-apertures (this
   # is computed from the pyramid WFS input during the calibration step).
   # The full list of parameters can be found in the init method of the PyrSlopec class.
   slopec:
     class:             'PyrSlopec'  
     pupdata_object:    'scao_pupdata'            # tag of the pyramid WFS pupils
     inputs:
       in_pixels:        'detector.out_pixels'

   # The modal reconstruction block reconstructs the wavefront slopes into modal coefficients.
   # It uses a reconstruction matrix that is computed during the calibration phase.
   rec:
     class:              'Modalrec'  
     recmat_object:      'scao_pyr_rec'         # reconstruction matrix tag
     inputs:
       in_slopes:        'slopec.out_slopes'
     outputs:  ['out_modes'] 

   # The control block computes the control commands based on the differential modal coefficients.
   # The modal coefficients are differential because it operates in closed loop.
   # The full list of parameters can be found in the init method of the Integrator class.
   control:
     class:             'Integrator'
     delay:             2                      # Total temporal delay in time steps
     int_gain:          [0.5]                  # Integrator gain (for 'INT' control)
     n_modes:           [54]                   # This means we use 54 modes with the same gain
     inputs:
       delta_comm: 'rec.out_modes'

   # The DM block simulates the deformable mirror.
   # In this case, it uses Zernike modes directly without generating a modal basis.
   # The Zernike modes are generated on the fly.
   # It can also use influence functions and modes-to-command matrix stored in files.
   # The DM height is set to 0, meaning it is at the pupil plane.
   # The full list of parameters can be found in the init method of the DM class.
   dm:
     class:             'DM'
     simul_params_ref:  'main'
     type_str:          'zernike'              # modes type
     nmodes:            54                     # number of modes
     obsratio:          0.0                    # obstruction dimension ratio w.r.t. diameter
     height:            0                      # DM height [m]
     inputs:
       in_command: 'control.out_comm'

   # The PSF block computes the point spread function (PSF) based on the electric field
   # at the pupil plane after the DM. It uses the wavelength and the padding coefficient
   # to compute the PSF. The start_time is the time at which the PSF integration starts.
   psf:
     class:             'PSF'
     simul_params_ref:  'main'
     wavelengthInNm:    1650                 # [nm] Imaging wavelength
     nd:                8                    # padding coefficient for PSF computation
     start_time:        0.05                # PSF integration start time
     inputs:
       in_ef:  'prop.out_on_axis_source_ef'
     outputs:  ['out_psf', 'out_sr']

   # Data store for saving the simulation results.
   # The data will be stored in a directory named with a timestamp (TN) located in 'output'.
   # In this case it saves the residual electric field on-axis source.
   # Other results can be added as needed.
   data_store:
     class:             'DataStore'
     store_dir:         './output'             # Data result directory: store_dir'/TN/'
     inputs:
       input_list: ['res_ef-prop.out_on_axis_source_ef']

If the full-rate data volume is too large, ``DataStore`` can downsample what is
written to disk for selected inputs:

.. code-block:: yaml

   data_store:
     class:             'DataStore'
     store_dir:         './output'
     downsample_factor_by_input:
       res_ef:          50
       sr:              5
     inputs:
       input_list: ['res_ef-prop.out_on_axis_source_ef',
                    'sr-psf.out_sr']

In this example, ``res_ef`` is stored every 50 received samples and ``sr``
every 5. The keys in ``downsample_factor_by_input`` are the aliases before the dash
in ``input_list``. If you want a single cadence for all inputs, use
``downsample_factor`` instead. The two options cannot be combined.

Part 2: Calibration
-------------------

Not all elements of an adaptive optics simulation are fully defined a priori.
Some key components—such as the geometry of valid subapertures for the wavefront sensor or the reconstruction matrix for modal control—depend on the specific configuration and must be determined through a **calibration process**.

Calibration in SPECULA consists of running dedicated, simplified simulations whose purpose is to "probe" the system and extract the necessary information for later use in the main (closed-loop) simulation.
These calibration runs typically:

- Use a reduced or modified version of the full simulation (for example, with only the pupil mask and no atmospheric or DM layers).
- Apply known inputs (such as push-pull commands or uniform illumination) to measure the response of specific components.
- Save the results (such as the valid pupil geometry or the interaction/reconstruction matrices) to files.

These calibration files are then referenced in the main YAML configuration and used by the corresponding objects during the closed-loop simulation.

This approach ensures that the simulation accurately reflects the real system's behavior.
Note that if you change the pupil geometry, WFS parameters, or DM configuration, you simply repeat the relevant calibration steps before running the full simulation.

In this case we need to calibrate two components:

1. **Pyramid WFS pupil geometry**: This defines the valid subapertures for the Pyramid WFS based on the pupil geometry.
2. **Reconstruction matrix**: This defines how the wavefront slopes measured by the Pyramid WFS are converted into modal coefficients for control. The reconstruction matrix is the inverse of the interaction matrix, which is computed by applying a known push-pull signal to the DM and measuring the resulting slopes.

Step 1: Calibrate the Pyramid WFS Pupil Geometry
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Create a YAML file, for example ``params_scao_pyr_test_pupdata.yml``:

.. code-block:: yaml

   # this is the block that generates the pupil geometry file used by the Pyramid WFS
   # as a function of a single frame of the Pyramid WFS
   pyr_pupdata:
     class: 'PyrPupdataCalibrator'
     thr1: 0.1
     thr2: 0.25
     display_debug: True
     output_tag:        'scao_pupdata'
     inputs:
       in_i: 'pyramid.out_i'

   # Override propagation parameters leaving only the pupilstop layer and removing
   # the atmo layers and the DM layer.
   # This is necessary to ensure the Pyramid WFS pupil geometry is correctly computed
   # considering the pupil mask.
   prop_override:
     inputs:
       common_layer_list: ['pupilstop']

   # Override main simulation parameters reducing the total time to 0.001s,
   # which is sufficient for getting a single frame of the Pyramid WFS.
   main_override:
     total_time:  0.001

   # Override Pyramid WFS parameters to set a very large modulation radius
   # This is necessary to reduce the diffraction effects and get uniform illumination
   pyramid_override:
     mod_amp:           10.0                     # Modulation radius (in lambda/D units)

   # Remove unnecessary components to avoid issues and speed up the calibration
   remove: ['atmo', 'slopec', 'rec','control','dm', 'data_store']

Run this calibration step with:

.. code-block:: bash

   specula params_scao_pyr_basic.yml params_scao_pyr_test_pupdata.yml

This will generate the ``scao_pupdata.fits`` file, which contains the valid pupil geometry for the Pyramid WFS.

Step 2: Calibrate the Reconstruction Matrix
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Create a YAML file, for example ``params_scao_pyr_test_rec.yml``:

.. code-block:: yaml

   # this block generates a positive and negative push-pull signal
   # to command the DM and compute the interaction matrix
   pushpull:
     class:     'PushPullGenerator'
     nmodes:    54
     vect_amplitude: [50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0,
                      50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0,
                      50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0,
                      50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0,
                      50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0,
                      50.0, 50.0, 50.0, 50.0] # Push-pull amplitude in nm (written by hand,
                      # it can be stored in a file)
     outputs:   ['output']

   # This block computes the interaction matrix using the Pyramid WFS
   # It takes the slopes from the Pyramid WFS and the push-pull commands as inputs
   im_calibrator:
     class:     'ImCalibrator'
     nmodes:    54
     im_tag:    'scao_pyr_im'
     data_dir:  './calib/im'
     overwrite: True
     inputs:
       in_slopes:   'slopec.out_slopes'
       in_commands: 'pushpull.output'
     outputs: ['out_intmat']

   # This block computes the reconstruction matrix from the interaction matrix
   # doing a pseudo-inverse operation
   rec_calibrator:
     class:     'RecCalibrator'
     nmodes:    54
     rec_tag:   'scao_pyr_rec'
     data_dir:  './calib/rec'
     overwrite: True
     inputs:
       in_intmat:   'im_calibrator.out_intmat'

   # Override main simulation parameters to set the total time
   # 54 times 2 (push+pull) times 0.001s
   main_override:
     total_time: 0.108                        # 54 modes × 2 (push+pull) × 0.001s

   # Override propagation parameters to include the pupilstop and DM layers
   # removing the atmo layers
   prop_override:
     inputs:
       common_layer_list: ['pupilstop', 'dm.out_layer']

   # Override the DM parameters to set the sign (by default it is -1)
   # and change the input command to the push-pull output
   dm_override:
     sign: 1
     inputs:
       in_command: 'pushpull.output'

   # Override the detector parameters to disable photon and readout noise
   detector_override:
     photon_noise:   False
     readout_noise:  False

   # Remove unnecessary components to avoid issues and speed up the calibration
   remove: ['atmo', 'rec','control']

Run this calibration step with:

.. code-block:: bash

   specula params_scao_pyr_basic.yml params_scao_pyr_test_rec.yml

This will generate the interaction matrix (``scao_pyr_im.fits``) and the reconstruction matrix (``scao_recmat.fits``).

Summary of Calibration Steps
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

1. **Pyramid Pupil Geometry:**  
   - Generates ``scao_pupdata.fits`` for the Pyramid WFS geometry.
2. **Reconstruction Matrix:**  
   - Generates ``scao_recmat.fits`` for modal reconstruction.

Both files are then referenced in your main YAML for closed-loop operation:

.. code-block:: yaml

   slopec:
     class:             'PyrSlopec'  
     pupdata_object:    'scao_pupdata'
     inputs:
       in_pixels:        'detector.out_pixels'

   rec:
     class:              'Modalrec'  
     recmat_object:      'scao_recmat'
     inputs:
       in_slopes:        'slopec.out_slopes'

---

Now your SCAO simulation is properly calibrated and ready for closed-loop runs!


Part 3: Running the Simulation
------------------------------

To run the simulation, use:

.. code-block:: bash

   specula params_scao_pyr_basic.yml

The simulation will run and save results in the output directory (``./output/``), including the residual electric field on-axis source as specified in the YAML file.

**Congratulations!**  
You have now configured and run a minimal SCAO simulation with a Pyramid WFS, Zernike DM, and a single atmospheric layer.

.. seealso::

   - :ref:`scao_tutorial` for a full-featured SCAO simulation, calibration and analysis workflow
   - :ref:`field_analyser_tutorial` for post-processing and analysis