Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5

ULySS: Tutorial on the creation of a new model component

The main strength of ULySS is its flexibility and extensibility. In particular, beside the model components proposed in the package, it is relatively easy for the programmer to provide new models.
ULySS is fitting a spectrum against a non-linear model. In fact this model is a linear combination of (any number of) arbitrary components. The purpose of this tutorial is to build a new component...
Basically, a component of the model fitted by ULySS is a function of any number of parameters returning a spectrum. ULySS will find the optimal set of values for these parameters in order to get the best match between the model and the spectrum to analyze.
In this tutorial we will construct a very simple component: A power-law spectrum that can be used to fit the nebular continuum present in the spectrum of an active galaxy.

To follow this tutorial it is advisable to study before the general tutorial, and possibly the other tutorials. It requires some knowledge of GDL/IDL (structures, pointers).

1. Create a model component for a power-law spectrum

In general a model component is described in a cmp structure described in the documentation of the function ULY_FIT. This structure contains some meta-data describing the component, but contains also the name and arguments of the functions to initialize and to evaluate the component (respectively the tags init_func, init_data and eval_fun, eval_data).
The definition of a component (ie. the creation of the cmp struture is made before calling the fitting procedure. For example, the function ULY_STAR defines a simple spectrum component (no free parameter).
The initialization of a component is performed when the sampling and wavelength range is known. It consists in doing everything that can be made only once, like reading a grid of models. init_fun is called by ULY_FIT_INIT with the argument init_data, which can be an arbitrary structure used only by eval_fun. The evaluation function produces eval_data stored in the component structure.
The evaluation is made at each step of the Levenberg-Marquart minimization and eval_fun is called with the argument eval_data and the array of parameters.

In the present case, the model-spectrum will be something like exp(findgen(npix)*step*tau), where step and npix are the logarithmic step and the number of pixels and tau the exponent of the power law. (note that the models are logarithmically sampled in order to fit the line-of-sight velocity distribution). There is a single non-linear parameter to fit, tau, and there is no need for an initialization, i.e. init_fun can be left as a null pointer.
The power-law component can be defined as: 
GDL> ; power is the component structure describing the power-low component
GDL> power = {name:'power', $
GDL>          descr:'Power-law component -- demonstration', $
GDL>          init_fun:'', $
GDL>          init_data:ptr_new(), $
GDL>          eval_fun:'tuto_power', $
GDL>          eval_data:ptr_new(/ALLO), $ ; will be filled by ULY_FIT_INIT
GDL>          para:ptr_new(/ALLO), $
GDL>          start:0d, $              ; this will be filled by ULY_FIT_INIT
GDL>          step:0d, $               ; this will be filled by ULY_FIT_INIT
GDL>          npix: 0l, $              ; this will be filled by ULY_FIT_INIT
GDL>          sampling:-1s, $          ; this will be filled by ULY_FIT_INIT
GDL>          mask:ptr_new(), $
GDL>          weight:0d, $
GDL>          e_weight:0d, $
GDL>          l_weight:0d, $
GDL>          lim_weig:[!VALUES.D_NaN,!VALUES.D_NaN] $
GDL>         }
GDL> 
GDL> ; We have also to define the para structure, describing the parameter tau:
GDL> *power.para = {name:'tau', unit:'', guess:ptr_new(0D), step:1D-2, $
GDL>                        limits:[-10d,10d], limited:[1,1], fixed:0S, $
GDL>                        value:0D, error:0D, dispf:''}
GDL> 
GDL> ; Since init_func is empty, ULY_FIT_INIT will fill the WCS information
GDL> ; with the values corresponding to the spectrum to fit, and will
GDL> ; pass the cmp structure as eval_data.
And the evaluation function as: 
function tuto_power, eval_data, pars
; pars is the array of parameters, actually only one value: the exponent
step = eval_data.step
npix = eval_data.npix

; simply return the array:
return, exp(findgen(npix)*step*pars[0])

end

2. Test the fit with a power-law spectrum

Lets first fit a galaxy spectrum with this component: 
GDL> ; Choose a galaxy spectrum to test this model, and try the fit
GDL> galaxy = uly_root+'/data/VazMiles_z-0.40t07.94.fits'
GDL> ulyss, galaxy, power, KMOMENT=0, MD=0, /PLOT

--------------------------------------------------------------------
INPUT PARAMETERS
--------------------------------------------------------------------
The fits file to be analyze is       ./data/VazMiles_z-0.40t07.94.fits
Name of the output file              output.res
Degree of multiplicative polynomial         0
No additive polynomial  
Component1 (power) Power-law component -- demonstration
  Guess for tau: 0.0000000 []
--------------------------------------------------------------------
SDSS style
Assume that the sampling is linear in wavelength
The error spectrum does not contain finite values
... ignore it (ie. do as if no errors were given)
--------------------------------------------------------------------
PARAMETERS PASSED TO ULY_FIT
--------------------------------------------------------------------
Wavelength range used                 :       3540.3541       7410.6865
 [Angstrom]
Sampling in log wavelength            :       51.501300 [km/s]
Number of independent pixels in signal:        4300
Number of pixels fitted               :        4300
DOF factor                            :      1.00000
--------------------------------------------------------------------
number of model evaluations:          16
time=       1.6582990
Number of pixels used for the fit        4296
----------------------------
chi^2          :       328.40524
----------------------------
Weight         :        193.29153 +/-      0.029005502 [data_unit/cmp_unit]
tau            :      0.80979726 +/- 0.00072975146 
----------------------------
As expected, the fit is poor ... but the program is working!
We can now simulate a composite spectrum, stellar population + power-law continuum, and analyse it with a composite model: 
GDL> g = uly_spect_read(uly_root+'/data/VazMiles_z-0.40t07.94.fits')
GDL> g = uly_spect_logrebin(g, /OVERWRITE)
GDL> tau = -3d
GDL> p = exp(findgen(n_elements(*g.data))*g.step*tau)
GDL> *g.data += 0.5 * mean(*g.data)/mean(p) * p 
GDL> 
GDL> window, 0
GDL> uly_spect_plot, galaxy, LINECOLOR='Black'
GDL> uly_spect_plot, g, /OVERPLOT
GDL> 
GDL> window, 1
GDL> cmp = [power, uly_ssp(MODEL_FILE=uly_root+'/models/Vaz_Miles.fits')]
GDL> ulyss, g, cmp, KMOMENT=0, MD=0, /PLOT

SDSS style
Assume that the sampling is linear in wavelength
The error spectrum does not contain finite values
... ignore it (ie. do as if no errors were given)
--------------------------------------------------------------------
INPUT PARAMETERS
--------------------------------------------------------------------
Name of the output file              output.res
Degree of multiplicative polynomial         0
No additive polynomial  
Component1 (power) Power-law component -- demonstration
  Guess for tau: 0.0000000 []
Component2 (cmp1) model:./models/Vaz_Miles.fits 
  Guess for age: 8000.0000 [Myr], Fe/H: -0.40000000 [dex]
--------------------------------------------------------------------
Read the model ./models/Vaz_Miles.fits
read the file ./models/Vaz_Miles.fits
model reading time=      0.66786098
model deriving time=      0.20855904
--------------------------------------------------------------------
PARAMETERS PASSED TO ULY_FIT
--------------------------------------------------------------------
Wavelength range used                 :       3540.3541       7410.6865
 [Angstrom]
Sampling in log wavelength            :       51.501300 [km/s]
Number of independent pixels in signal:        4300
Number of pixels fitted               :        4300
DOF factor                            :      1.00000
--------------------------------------------------------------------
number of model evaluations:          66
time=       1.2193091
Number of pixels used for the fit        4296
----------------------------
chi^2          :   1.2476216e-22
----------------------------
Light fraction :        33.333333 +/-      0.018466303 %
Weight         :        328.21149 +/-       0.18182558 [data_unit/cmp_unit]
tau            :      -3.0000000 +/- 0.017293929 
----------------------------
Light fraction :        66.666667 +/-      0.020968968 %
Weight         :       1.00000000 +/-    0.00031453452 [data_unit/cmp_unit]
age            :       7940.0005 +/- 80.762597 Myr
Fe/H           :     -0.40000001 +/- 0.0054215726 dex
----------------------------
This new component is working ... before using it in astronomical applications it will have to be tested carefully (study the degeneracies and the local minima).
Contact: ulyss at obs.univ-lyon1.fr Last modified: Thu Jul 31 13:57:09 2008.