.. SPDX-FileCopyrightText: 1992-2026 NWO-I/SRON Space Research Organisation Netherlands
..
.. SPDX-License-Identifier: CC-BY-4.0

.. _spexapi_fit:

Optimization & Simulation functions
===================================

Fitting spectra
---------------

.. _spexapi_eval:

Evaluate the model and convolve with response
"""""""""""""""""""""""""""""""""""""""""""""

The smallest step in a fitting process is calculating the model,
fold the model spectrum through the response matrix, and
calculate the likelihood compared to a loaded dataset.

The eval procedure calculates the model and folds it
through the response matrix::

    call spex%eval()

The interface looks like this::

      module subroutine spex_fit_eval(this)
        class(sapi), intent(inout)  :: this
      end subroutine

If you need the fit statistics and degrees of freedom, they
can be returned using the ``eval_stat`` procedure::

    call spex%eval_stat(stat,dof)

The interface looks like this::

      module subroutine spex_fit_eval_stat(this, stat, dof)
        class(sapi),  intent(inout)  :: this
        real, intent(out)            :: stat   !! Statistics value
        integer, intent(out)         :: dof    !! Degrees of freedom
      end subroutine

.. _spexapi_fit:

Fitting
"""""""

Spectral fitting is done with the ``fit`` procedure::

    call spex%fit(ier)

    or

    call spex%fit(ier,iter=10)

The interface looks like this::

      module subroutine spex_fit_fit(this, ier, iter)
        class(sapi),  intent(inout)      :: this
        integer, intent(inout)           :: ier   !!Error status
        integer, intent(inout), optional :: iter  !!Number of iterations
      end subroutine

One can tweak the fits with the commands below. The fit method can be set with::

    call spex%fit_method('classical',ier)

Interface::

      module subroutine spex_fit_set_method(this, imethod, ier)
        class(sapi),  intent(inout)    :: this
        character*(*),  intent(in)     :: imethod      !! Method fitting
        integer, intent(out)           :: ier          !! Error status
      end subroutine

The terminal output for each iteration can be set with::

    call spex%fit_print(1,ier)

Interface::

      module subroutine spex_fit_set_print(this, i, ier)
        class(sapi),  intent(inout)   :: this
        integer, intent(in)           :: i       !!printing the intermediate step (yes if i=1)
        integer, intent(out)          :: ier     !!error status
      end subroutine

The fit statistics can be set using::

    call spex%fit_stat('cstat',ier)

Interface::

      module subroutine spex_fit_set_statistic(this, stat, ier)
        class(sapi),  intent(inout)    :: this
        character*(*), intent(in)      :: stat     !! Fit statistics (chi2, cstat, or wstat)
        integer, intent(out)           :: ier      !! error status
      end subroutine

And setting the statistics for each instrument separately::

    call spex%fit_stat_inst('cstat',1,1,ier)

Interface::

      module subroutine spex_fit_set_statistic_inst(this, stat, inst, ireg, ier)
        class(sapi),  intent(inout)    :: this
        character*(*), intent(in)      :: stat      !! Fit statistics (chi2, cstat, or wstat)
        integer, intent(in)            :: inst      !! Instrument number
        integer, intent(in)            :: ireg      !! Region number
        integer, intent(out)           :: ier       !! error status
      end subroutine

There is also a deprecated option to set the weigths for the fit::

    call spex%fit_weight('model',ier)

Interface::

      module subroutine spex_fit_set_weight(this, fweight, ier)
        class(sapi),  intent(inout)    :: this
        character*(*),  intent(in)     :: fweight  !!fit weight: 'data' or 'model'
        integer, intent(out)           :: ier      !!error status
      end subroutine


.. _spexapi_error:

Error calculation
-----------------

The error on a parameter can be calculated with::

    call spex%error(isect,icomp,pnam,pvalue,perr_low,perr_upp,perr_lc,perr_cmin,perr_pmin,ier)

Interface::

      module subroutine spex_fit_error(this, isect, icomp, pnam, pvalue, perr_low, perr_upp, &
                                       perr_lc, perr_cmin, perr_pmin, ier)
        class(sapi),  intent(inout)     :: this
        integer, intent(in)             :: isect    !! The instrument number
        integer, intent(in)             :: icomp    !! The component number
        character*(*), intent(inout)    :: pnam     !! Parameter name
        real(dp), intent(out)           :: pvalue   !! Parameter value
        real(dp), intent(out)           :: perr_low !! Lower error
        real(dp), intent(out)           :: perr_upp !! Upper error
        logical, intent(out)            :: perr_lc  !! Lower chi/c statistics found (True/False)
        real(dp), intent(out)           :: perr_cmin!! Lowest chi/c statistics found
        real(dp), intent(out)           :: perr_pmin!! Parameter value for new minimum
        integer, intent(inout)          :: ier      !! status error
      end subroutine

Example::

    call spex%error(1,1,'norm',pvalue,perr_low,perr_upp,perr_lc,perr_cmin,perr_pmin,ier)

One can also change the target delta-chi::

    call spex%error_dchi(dchi)

Interface::

      module subroutine spex_fit_set_error_dchi(this,dchi)
        class(sapi),  intent(inout)  :: this
        real, intent(in)             :: dchi !! Delta chisq of the error to search for
      end subroutine

Example::

    call spex%error_dchi(1.0)

Setting the error start value can be done with::

    call spex%error_start(start)

Interface::

      module subroutine spex_fit_set_error_start(this,start)
        class(sapi),  intent(inout) :: this
        real, intent(in)            :: start !! Starting value of error calculation
      end subroutine

Example::

    call spex%error_start(1E+5)

.. _spexapi_step:

Step
----

Similar to the step command in SPEX, the step process is set up with multiple
procedures::

    call spex%step_axis(iaxis,isect,icomp,npar,lim1,lim2,nstep,log)

The above sets up one stepping axis for a parameter, including the limits and stepsize.

The interface::

      module subroutine spex_fit_step_axis(this, iaxis, isect, icomp, npar, lim1, lim2, nstep, log)
        class(sapi),  intent(inout) :: this
        integer, intent(in)     :: iaxis    !! Axis number to set
        integer, intent(in)     :: isect    !! Sector number
        integer, intent(in)     :: icomp    !! Component number
        character*4, intent(inout) :: npar     !! Parameter name
        real, intent(in)        :: lim1     !! Lower limit of axis
        real, intent(in)        :: lim2     !! Upper limit of axis
        integer, intent(in)     :: nstep    !! Number of steps
        logical, intent(in)     :: log      !! Axis logarithmic (yes=.true., no=.false.)
      end subroutine

Check if the axes are properly defined::

    call spex%step_check_axis(isect,icomp,npar,lim1,lim2,nstep,log,ier)

If ``ier`` is non-zero, the axes were not set up properly.

The interface::

      module subroutine spex_fit_step_check_axis(this, isect, icomp, npar, lim1, lim2, nstep, log, ier)
        class(sapi),  intent(inout) :: this
        integer, intent(in)     :: isect    !! Sector number
        integer, intent(in)     :: icomp    !! Component number
        character*4, intent(inout) :: npar     !! Parameter name
        real, intent(in)        :: lim1     !! Lower limit of axis
        real, intent(in)        :: lim2     !! Upper limit of axis
        integer, intent(in)     :: nstep    !! Number of steps
        logical, intent(in)     :: log      !! Axis logarithmic (yes=.true., no=.false.)
        integer, intent(inout)  :: ier
      end subroutine

Example::

    call spex%step_check_axis(1,1,'norm',)

Set the number of axes for the step session::

    call spex%step_dim(ndim,ier)

Interface::

      module subroutine spex_fit_step_dim(this, ndim, ier)
        class(sapi),  intent(inout)   :: this
        integer, intent(in)           :: ndim      !! Step dimension
        integer, intent(out)          :: ier
      end subroutine

Setting an output file for step is done with::

    call spex%step_file('stepname')

Interface::

      module subroutine spex_fit_step_file(this, stpname)
        class(sapi),  intent(inout) :: this
        character*(*), intent(in)   :: stpname
      end subroutine

The final step is started with::

    call spex%step(ier)

Interfaces::

      module subroutine spex_fit_step(this,ier,pstep)
        class(sapi),  intent(inout)   :: this
        integer, intent(inout)        :: ier       !error status
        type(step), optional, pointer :: pstep
      end subroutine

      module subroutine spex_fit_step_check(this, ier, pstep)
        class(sapi),  intent(inout)   :: this
        integer, intent(inout)        :: ier
        type(step), pointer, optional :: pstep
      end subroutine


.. _spexapi_sim:

Simulating spectra
------------------

Following the simulate command in SPEX, the Fortran API has
procedures to set the different simulation properties if needed.

Simulating noise on the background or source can be switched on
and off::

    call spex%sim_noise(.false.)
    call spex%sim_bnoise(.false.)

The interfaces for these procedures::

      module subroutine spex_fit_sim_noise(this,noise)
        class(sapi),  intent(inout)   :: this
        logical, intent(in)          :: noise     !! Set source noise for simulation (true/false)
      end subroutine

      module subroutine spex_fit_sim_bnoise(this,bnoise)
        class(sapi),  intent(inout)   :: this
        logical, intent(in)          :: bnoise   !! Set background noise for simulation (true/false)
      end subroutine

The random seed can be set to a value or randomized with::

    call spex%sim_seed(.true.,0)          ! Set random seed to a random number based on computer clock
    call spex%sim_seed(.false.,12345678)  ! Set random seed to a particular number

The interface::

      module subroutine spex_fit_sim_seed(this,rand,seed)
        class(sapi),  intent(inout)   :: this
        logical, intent(in)           :: rand   !! Should the seed be randomly set? True/False
        integer, intent(in)           :: seed   !! Set seed manually if rand is False
      end subroutine

Additional systematic errors can be added (not recommended)::

    call spex%sim_syserr(sys_src,sys_bkg)

Interface::

      module subroutine spex_fit_sim_syserr(this,sys_src,sys_bkg)
        class(sapi),  intent(inout)    :: this
        real, intent(in)               :: sys_src  !! Systematic error on the source spectrum
        real, intent(in)               :: sys_bkg  !! Systematic error on the background spectrum
      end subroutine

Select the instruments that need to be simulated::

    call spex%sim_inst(in1,in2,ic1,ic2,ier)

The interface::

      module subroutine spex_fit_sim_inst(this,in1,in2,ic1,ic2,ier)
        class(sapi), intent(inout)    :: this
        integer, intent(in)           :: in1
        integer, intent(in)           :: in2
        integer, intent(in)           :: ic1
        integer, intent(in)           :: ic2
        integer, intent(inout)        :: ier
      end subroutine

And finally, the main simulate command::

    call spex%sim(exposure)

The interface::

      module subroutine spex_fit_simulate(this, exposure)
        class(sapi),  intent(inout)   :: this
        real(dp), intent(in)          :: exposure   !! Exposure time
      end subroutine

Example::

    call spex%sim(1E+6)