# document.write(pagetitle)

Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:
;                  ULYSS
; PURPOSE:
;                  Analyse a spectrum
; USAGE:
;                  ulyss,
;                  <spectrum>,<cmp> or MODEL_FILE=<model_file>
;                  [, POSITION=<position>]
;                  [, ERR_SP=<err_sp>][, SNR=<snr>]
;                  [, VELSCALE]
;                  [, KMOMENT=<kmoment>][, SG=<sg>][, DG=<dg>][, DL=<dl>]
;                  [, LMIN=<lmin>][, LMAX=<lmax>]
;                  [, NSIMUL=<nsimul>]
;                  [, KFIX=<kfix>][, KPEN=<kpen>]
;                  [, /CLEAN]
;                  [, /QUIET]
;                  [, SOLUTION=<solution>]
;                  [, /ALL_SOLUTIONS]
;                  [, FILE_OUT=<result_file>]
;                  [, MODECVG=<modecvg)]
;                  [, /PLOT]
;
; DESCRIPTION:
;      Main procedure of the ULySS package.
;      Reads the observed spectrum, the models, and makes the fit (with uly_fit).
;
;      ULYSS fits a spectrum with a linear combination of non-linear components
;      convolved with a given line-of-sight velocity distribution (LOSVD)
;      and multiplied by a polynomial continuum. The linear or non-linear
;      parameters may be bounded (for example force a non-negative
;      combination of components).
;
;      The fit is a Levenberg-Marquart local minimization for the non-linear
;      parameters and a bounded-values least square for the linear ones.
;      Any parameter can be bounded of fixed. The algorithm is described
;      in ULY_FIT.
;
;      Line-of-sight velocity distribution (LOSVD):
;
;      Velocity dispersion:
;      The velocity dispersion computed by the program is
;      sigma_ulyss = sqrt(sigma_obs^2 - sigma_model^2). Where
;      sigma_obs = sqrt(sigma_physical^2 + sigma_instrumental^2) and
;      sigma_model is the dispersion in the model. For ELODIE-based models,
;      sigma_model is 13 km/s; For MILES-based models, it is ~58km/s.
;      sigma_instrumental is the instrumental broadening of the observation.
;      If the relative line-spread function of the spectrograph is injected in
;      the model, then sigma_ulyss = sigma_phys.
;
;      Logarithmic sampling:
;      The line of sight velocity distribution is measured as a broadening of
;      the spectral features due to Doppler shift. The wavelength shift is
;      related to the redshift z by: (l1-l0)/l0 = z, where l1 and l0 are
;      respectively the shifted and restframe wavelengths. (If we neglect the
;      relativistic term, z=v/c, where v is the velocity shift and c the speed
;      of light). As the wavelength shift depends on the wavelength, when
;      sampled linearly in wavelength, a spectrum is not simply 'translated'.
;      However, with a logarithmic sampling the shift becomes a translation,
;      and the effect of the LOSVD can be written as a convolution.
;
;      Note also to take care of the composition of the partial shift:
;      If a shift z1 is applied in ULY_SPECT_READ, and if ULYSS finds
;      a 'residual' shift z2, the total redshift is:
;      1 + z = (1 + z1)(1 + z2), when z1xz2 is small it is close to 1+z1+z2.
;
;      Relativistic correction to compute the radial velocity:
;
;      In special relativity, the relation between z and v is:
;      1 + z = sqrt(c+v/c-v) = (1 + v/c) / sqrt(1-(v/c)^2).
;      (The last term is the Lorentz factor, gamma: 1+z = gamma (1+v/c)).
;      At low redshift, the Newtonian approximation is: z=v/c
;
;      This formula can be inverted into:
;      v/c = ((1+z)^2 - 1) / ((1+z)^2 + 1)
;      (z=1 corresponds to v = 0.6 c = 180000 km/s)
;
;      ULySS determines cz, and no relativistic correction is applied.
;      The transformation above must be used to determine v from cz.
;
;      In many cases computing cz is sufficient. The databases, LEDA or NED
;      give cz that is somehow inproperly sometime called a 'velocity'.
;
;      Computation of the velocity dispersion:
;
;      ULySS determines the cz dispersion, and in general this shall be
;      converted into a *velocity* dispersion to carry on dynamics (for
;      example to apply the Virial Theorem).
;      Fortunately, sigma_cz and sigma are usually very similar.
;
;      To be rigorous, two effects shall be considered (i) the effect of the
;      composition of redshifts, and (ii) the relativistic correction.
;      Lets write the redshift of a star with v=sigma, in the galaxy's restframe
;      zs = (1 + sigma/c) / sqrt(1-(sigma/c)^2) - 1
;      And the residual redshift of the barycenter of the galaxy (after
;      de-redshifting the spectrum, for example using ULY_SPECT_READ): zb.
;      The cz dispersion is: sigma_cz = c (1+zb) zs
;      The two corrective terms are usually negligible at any redshift.
;      The first order correction pointed by Harrison (1974 ApJ, 191, L51)
;      for clusters of galaxies does not apply here.
;
; ARGUMENTS:
;   <spectrum>     Input
;                  File name, or 'spect' structure containing 1d spectrum
;                  to be analysed with uly_fit.
;                  A 'spect' structure is returned by, e.g. ULY_SPECT_READ.
;   <cmp>          Input
;                  Array of model's components, see ULY_FIT.
;                  Prior calls to ULY_STAR, ULY_TGM, ULY_SSP and/or ULY_LINE
;                  can define this array and set the guesses and constraints
;                  on the free parameters.
;                  MODEL_FILE can be given instead.
;
; KEYWORDS:
;
;   The following keywords are handled by ULYSS itself, or shared
;
;   MODEL_FILE     Input, filename
;                  When <cmp> is not provided, this keyword may give the
;                  name of a FITS file containing the model to fit.
;                  <cmp> and MODEL_FILE are mutually exclusive.
;                  Note that using MODEL_FILE is less flexible that using
;                  the component-definition functions. See ULY_SSP, ULY_TGM or
;                  other component-definition functions for more capabilities.
;
;   SNR            Input, float
;                  Mean signal to noise ratio the of analyzed spectrum.
;                  This parameter is used to derive the errors if the error
;                  spectrum is not attached to the input spectrum and if ERR_SP
;                  is not given.  It will generate a constant error spectrum.
;                  SNR is ignored if an error spectrum is available.
;
;   NSIMUL         Input, integer
;                  To make Monte-Carlo simulations, set with this keyword
;                  the number of simulations. These simulations are
;                  noise. To take into account the correlation of the noise
;                  introduced along the data processing (for example,
;                  resulting of an oversampling), the noise is generated
;                  on a vector of length npix/dof_factor, and then rebinned
;                  to npix. Where npix is the actual number of pixels, and
;                  dof_factor characterizes the correlation of the noise.
;                  (dof_factor is included in the spect structure,
;                  ULY_SPECT_LOGREBIN and the other smoothing or resampling
;                  routines modify it consistently).
;
;   /QUIET
;                  Set this keyword to supress the printing of information and
;                  results.
;
;   FILE_OUT
;                  The names of the result file are constructed by appending
;                  '.res' and '.fits' to this variable.
;                  The '.res' ASCII file contains the values of the parameters
;                  and their uncertainties. If the file pre-exists, new records
;                  are appended. This file can be used by ULY_SOLUT_TPRINT,
;                  ULY_SOLUT_TPLOT, ...
;                  The FITS file contains the spectrum, the bestfit, the
;                  polynomials and the mask of good pixels. In can be plotted
;                  with ULY_SOLUT_SPLOT.
;                  If neither FILE_OUT or SOLUTION are specified, output files
;                  with prefix 'output' are created
;
;   SOLUTION
;                  Output structure containing all fitted parameters and their
;                  respective errors. See ULY_FIT for details.
;
;   /ALL_SOLUTIONS
;                  When <cmp> specifies a grid of guesses (for global
;                  minimization), this keyword tells to return all the local
;                  solutions. By default, only the best solution is returned.
;
;   /PLOT
;                  Set this keyword to display the fit using ULY_SOLUT_SPLOT.
;
;   POSITION
;                  When a multidimensional dataset is to be analysed, like a
;                  long-slit spectrum, stacked spectra or cube, the keyword
;                  specifies the position of the 1D spectrum to analyse. It can
;                  be a scalar or an array of the same dimension as the dataset
;                  (not counting the spectral dimension).
;
;   VELSCALE       [km/s] default=conserve the number of pixels
;                  Sampling. Size of the pixel after the rebinning in logarithm
;                  of wavelength.
;
;   and are therefore relevant only if <spectrum> is a file name (if <spectrum>
;   is a spect structure, ULY_SPECT_READ is not called).
;   Check in the documentation of ULY_SPECT_READ for further information.
;
;   SG             [dimensionless] default=0
;                  Guess for the redshift z. This keyword is handled by
;                  uly_spect_read to shift the data. (the guess for the
;                  minimization is set to 0, and cannot by changed).
;                  This value must be quite precise: An error by more than 3 or
;                  5 times the velocity dispersion may prevent the fit to
;                  converge.
;
;   LMIN,LMAX      [Angstrom]
;                  Minimum and maximum wavelength. These parameters can
;                  be vectors to define several fitting intervals.
;
;   ERR_SP
;                  Name of a FITS file containing the error spectrum.
;                  (In some cases the error spectrum may be included in the
;                  same file as the signal, see ULY_SPECT_READ)
;
;   The following keywords are handled by ULY_FIT. Check its documentation
;   for further information.
;
;   KMOMENT
;                  Number of terms of the LOSVD.
;                  The terms are in the order: [cz, sigma, h3, h4, h5,
;                  h6]. By default, KMOMENT=2, i.e. a Gaussian LOSVD
;                  is fitted.
;
;   DG             [km/s]
;                  Guess or fixed value for the velocity dispersion.
;                  This guess is not very critical, and the default value
;                  of 1 pixel is generally satisfactory.
;
;   DL             [km/s]
;                  2-elements array giving the limits for the fitted velocity dispersion.
;
;   KFIX
;                  Array used to fix some of the parameters of the LOSVD,
;                  0 means that the parameter is free; 1 that it is fixed.
;                  The parameters are specified in the following order:
;                  [cz, sigma, h3, h4, h5, h6]. cz and sigma must be
;                  given in km/s.
;
;   KPEN
;                  Kinematics penalty parameter for the Gauss-Hermite
;                  coefficients. Default is no penalization (KPEN=0).
;                  This penalisation is described in Cappellari & Emmsellem
;                  2004, PASP 116, 138. The actual value should be chosen
;                  carefully (performing Monte-Carlo simulations), but values
;                  in the range 0.5 to 1.0 shall be a good starting guess.
;
;                  Default -1; Degree of additive polynomial,
;
;   MDEGREE
;                  Default 10; Degree of multiplicative polynomial
;
;   POLPEN
;                  Automatic selection of significant terms in the
;                  multiplicative polynomial. See ULY_FIT_LIN.
;
;   CLEAN
;                  Set this keyword to iteratively detect and clip the outiers.
;                  The algorithm is described in ULY_FIT.
;
;   MODECVG
;                  Keyword passed to ULY_FIT (see the documentation there)
;
; EXAMPLE:
;   1. Simple and quite limited examples:
;     Lets first select a spectrum file provided in the package for tests:
;     (the variable uly_root is set in uly_startup.pro to point to ULySS
;      root directory)
;           star = uly_root+'/data/cflib_114642.fits'
;
;     Analyse this file with default parameters (Determine the broadening
;     by comparing to the solar spectrum):
;           ulyss, star
;
;     Determine the atmospheric parameters of this star with the Elodie.3.2
;     interpolator
;           ulyss, star, MODEL_FILE=uly_root+'/models/elodie32_flux_tgm.fits'
;
;     Analyse a stellar population with a Pegase_HR/Elodie.3.1 grid of SSPs:
;           galaxy = uly_root+'/data/VazMiles_z-0.40t07.94.fits'
;           ulyss, galaxy, MODEL_FILE=uly_root+'/models/PHR_Elodie31.fits'
;
;     Analyse a spectrum previously read and convolved
;           galaxy = uly_root+'/data/VazMiles_z-0.40t07.94.fits'
;           spectrum = uly_spect_losvdconvol(spectrum, 0., 30., 0, 0, /OVER)
;           ulyss, spectrum, MODEL=uly_root+'/models/PHR_Elodie31.fits'
;
;     In real cases, the model should be defined before calling ULYSS,
;     in order to choose the guesses, add constraints on the parameters,
;     or define a composite model.
;
;     Analyse a stellar population with Pegase_HR/Elodie.3.1
;           galaxy = uly_root+'/data/VazMiles_z-0.40t07.94.fits'
;           cmp1 = uly_ssp(AG=[1000.], ZG=[0], AL=[200.,3000.])
;           cmp2 = uly_ssp(AG=[10000.], ZG=[0], AL=[3000., 12000.])
;           cmp = [cmp1, cmp2]
;           ulyss, galaxy, cmp, /PLOT
;
;     See the routines ULY_SSP, ULY_TGM, ULY_STAR, ...
;     to define model-components which can be assembled together as a vector.
;
; AUTHOR:
;                  Mina Koleva, 03/03/2006
;
; HISTORY:
;     2006/03/03    Created by Mina Koleva for stellar population fitting
;     2007/07/06    Change keywords,rewrite doc...
;     2008/05/27    Restructured and made to fit any model, Ph. Prugniel
;
;-
; CATEGORY:    ULY
; Copyright (C) 2008 Koleva M. , Prugniel Ph., Bouchard A., Wu Y.
;
;    This program is free software: you can redistribute it and/or modify
;    the Free Software Foundation, either version 3 of the License, or
;    (at your option) any later version.
;
;    This program is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;    GNU General Public License for more details.
;
;    You should have received a copy of the GNU General Public License
;    along with this program.  If not, see <http://www.gnu.org/licenses/>.
;
;------------------------------------------------------------------------------
; utility program to read cmp from a file
compile_opt idl2, hidden

if size(model_file, /TYPE) ne 7 then begin
message, 'Argument model_file must be a filename', /INFO
return, 0
endif
file = strtrim(model_file, 2)
if file_test(file) ne 1 then begin
file += '.fits'
if file_test(file) ne 1 then begin
message, 'Argument model_file must be a filename ('+model_file+')', /INFO
return, 0
endif
endif

if mess ne '' then begin
message, mess, /INFO
return, 0
endif

case strtrim(string(sxpar(h,'ULY_TYPE')),2) of
'SSP' : return, uly_ssp(MODEL_FILE=model_file)
'TGM' : return, uly_tgm(MODEL_FILE=model_file)
else: return, uly_star(model_file)
endcase

end

;------------------------------------------------------------------------------
pro ulyss, inspectr, cmp,                         \$
MODEL_FILE=model_file,                 \$
FILE_OUT=file_out,                     \$
POSITION=position,                     \$
ERR_SP=err_sp, SNR=snr,                \$
SG=shift_guess,                        \$
KMOMENT=kmoment, DG=sigma_guess,       \$
DL=sigma_limits,                       \$
VELSCALE=velscale,                     \$
MDEGREE=mdegree, POLPEN=polpen,        \$
LMIN=lmin, LMAX=lmax,                  \$
KFIX=kfix, KPEN=kpen,                  \$
NSIMUL=nsimul,                         \$
CLEAN=clean,                           \$
QUIET=quiet,                           \$
MODECVG=modecvg,                       \$
SOLUTION=solution,                     \$
ALL_SOLUTIONS=all_solutions,           \$
PLOT=plot

;%%%%%%%%%%%%% PART I - TESTS, INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
compile_opt idl2
on_error, 2

; check the input parameters and set the default values
;
if n_elements(inspectr) eq 0 then begin
file = ' '
read, file, prompt = 'give name of input file: '
endif else if size(inspectr,/TYPE) ne size('str',/TYPE) then begin
spectrum = inspectr
file = '<spectrum>'
;   check that neither SG, LMIN or LMAX are given
if n_elements(shift_guess) ne 0 then \$
message, 'ULYSS: When <spectrum> is a structure SG=sg should not be specified'
if n_elements(lmin) ne 0 or n_elements(lmax) ne 0 then \$
message, 'ULYSS: When <spectrum> is a structure LMIN and LMAX should not be specified'
if n_elements(err_sp) ne 0 then \$
message, 'ULYSS: When <spectrum> is a structure ERR_SP should not be specified'
endif else file = strtrim(inspectr,2)
if file ne '<spectrum>' then begin
testinp = file_test(file)
if testinp eq 0 then \$
message,'The file:'+file+' does not exist. Give correct name for the input file.'
endif

if n_elements(cmp) gt 0 then if n_elements(model_file) gt 0 then \$
message, 'Invalid arguments, <cmp> and MODEL_FILE are exclusive'

; Save the input arguments in 'command' for header of .res file
if file ne '<spectrum>' then command = 'ulyss, ''' + file + '''' \$
else command = 'ulyss, ' + file
if n_elements(file_out) ne 0 then command += ', FILE_OUT='''+file_out+''''
if n_elements(shift_guess) ne 0 then command += ', SG='+strtrim(shift_guess,2)
if n_elements(kmoment) ne 0 then command += ', KMOMENT='+strtrim(kmoment,2)
if n_elements(kpen) ne 0 then command += ', KPEN='+strtrim(kpen,2)
if n_elements(sigma_guess) ne 0 then command += ', DG='+strtrim(sigma_guess,2)
if n_elements(snr) ne 0 then command += ', SNR='+strtrim(snr,2)
if n_elements(velscale) ne 0 then command += ', VELSCALE='+strtrim(velscale,2)
if n_elements(mdegree) ne 0 then command += ', MDEGREE='+strtrim(mdegree,2)
if n_elements(lmin) gt 0 then begin
command += ', LMIN=['+strtrim(lmin[0],2)
for kk=1,n_elements(lmin)-1 do command += ','+strtrim(lmin[kk],2)
command += ']'
endif
if n_elements(lmax) gt 0 then begin
command += ', LMAX=['+strtrim(lmax[0],2)
for kk=1,n_elements(lmax)-1 do command += ','+strtrim(lmax[kk],2)
command += ']'
endif
if n_elements(kfix) gt 0 then begin
command += ', KFIX=['+strtrim(kfix[0],2)
for kk=1,n_elements(kfix)-1 do command += ','+strtrim(kfix[kk],2)
command += ']'
endif
if n_elements(nsimul) ne 0 then command += ', NSIMUL='+strtrim(nsimul,2)
if n_elements(err_sp) ne 0 then command += ', ERR_SP='+strtrim(err_sp,2)
if keyword_set(clean) ne 0 then command += ', /CLEAN'
if n_elements(polpen) ne 0 then command += ', POLPEN='+strtrim(polpen,2)
if keyword_set(quiet) ne 0 then command += ', /QUIET'

; Set default values
if n_elements(mdegree) ne 1 then mdegree = 10

if not arg_present(solution) and n_elements(file_out) eq 0 then file_out = 'output'

if n_elements(file_out) ne 0 then begin
file_res = file_out + '.res'
file_fits = file_out
endif

if n_elements(kmoment) eq 0 then kmoment = 2; moments of the gaussian to be fit

; kfix is an array of 0/1 specifying if a parameter of the LOSVD
; is fixed(1) or free(0)
if n_elements(kfix) gt 0 then begin
if n_elements(kfix) gt kmoment then \$
message, 'The number of elements of KFIX should not exceed '+\$
strtrim(string(kmoment),2)
endif

if n_elements(cmp) eq 0 then begin
if n_elements(model_file) gt 0 then begin
if size(cmp, /TYPE) ne 8 then message, 'Could not read model'
endif else begin
if not keyword_set(quiet) then \$
print, 'Use default model: Solar spectrum'
common uly_path, uly_root
cmp = uly_star(uly_root+'/models/sun.fits')
if size(cmp, /TYPE) ne 8 then \$
message, 'Could not define the model'
endelse
endif

if not keyword_set(quiet) then begin
print,'--------------------------------------------------------------------'
print, 'INPUT PARAMETERS'
print,'--------------------------------------------------------------------'
if file ne '<spectrum>' then \$
print, 'The fits file to be analyze is       ', file
if n_elements(file_res) gt 0 then \$
print, 'Name of the output file              ', file_res
print, 'Degree of multiplicative polynomial  ',mdegree
if adegree eq -1 then \$
print, 'No additive polynomial  ' \$
else \$

for n=0,n_elements(cmp)-1 do begin
print, 'Component',strtrim(n+1,2), ' (',cmp[n].name, ') ', cmp[n].descr
if n_elements(*cmp[n].para) gt 0 then print, '  Guess for ', FORM='(A,\$)'
outs = 0
for k=0,n_elements(*cmp[n].para)-1 do begin
if strlen((*cmp[n].para)[k].name) gt 0 then begin
if (*cmp[n].para)[k].dispf eq 'exp' then begin
guess = exp(*(*cmp[n].para)[k].guess)
endif else begin
guess = *(*cmp[n].para)[k].guess
endelse
if outs eq 1 then print, ', ', FORMAT='(A,\$)'
print, (*cmp[n].para)[k].name, ': ', \$
strjoin(strtrim(guess,2),' '), \$
' [', (*cmp[n].para)[k].unit,']', FORMAT='(6A,\$)'
outs = 1
endif
endfor
if n_elements(*cmp[n].para) gt 0 then print, '' ; end the line
endfor

print,'--------------------------------------------------------------------'
endif

c = 299792.458d                 ; [km/s] Speed of light

;%%%%%%%%%%%%%%%%%%%%  PART II  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%%%%%%%%%%% READ, RESAMPLE, CUT OBSERVATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
time0 = systime(1)

; SignalLog may be 2D (like a long-slit), in that case extract a scan
if n_elements(position) ne 0 then pos = position else pos = 0

if n_elements(spectrum) ne 0 then begin        ; if we pass a spectrum struct
SignalLog = uly_spect_extract(spectrum, ONED=pos, STATUS=status)
if status ne 0 then message, 'Could not extract 1D spectrum'
SignalLog = uly_spect_logrebin(SignalLog, velscale, /OVER)

endif else begin                ; if we pass a filename
SignalLog = uly_spect_read(file, lmin, lmax, VELSCALE=velscale, \$
ERR_SP=err_sp, SG=shift_guess, QUIET=quiet \$
)
SignalLog = uly_spect_extract(SignalLog, ONED=pos, /OVER, STATUS=status)
if status ne 0 then message, 'Could not extract 1D spectrum'
if SignalLog.sampling ne 1 then \$
SignalLog = uly_spect_logrebin(SignalLog, velscale, /OVER)
endelse

if not uly_spect_get(SignalLog, /VALID) then begin
message, 'Input spectrum is invalid', /CONT
return
endif
if n_elements(*SignalLog.goodpix) eq 1 then if *SignalLog.goodpix eq -1 then begin
message, 'There are no good pixels in input spectrum', /CONT
return
endif
if n_elements(nsimul) gt 0 then if n_elements(*SignalLog.err) eq 0 and n_elements(snr) eq 0 then begin
message, 'Cannot use NSIMUL because there is no noise spectrum', /CONT
return
endif

; make the good pixels list (if it does not exist) ... needed for the
; tests below
if n_elements(*SignalLog.goodpix) gt 0 then gp = *SignalLog.goodpix \$
else gp = lindgen((size(*SignalLog.data,/DIM))[0])

; compute the error if not read from the file (use SNR)
if n_elements(*SignalLog.err) eq 0 and n_elements(snr) ne 0 then begin
mean_error = mean((*SignalLog.data)[gp], /NAN) / snr
if not finite(mean_error) then \$
message, 'Cannot compute the mean of the signal!'
*SignalLog.err = (*SignalLog.data) * 0 + mean_error
endif

if n_elements(*SignalLog.err) gt 0 then begin
; Check if the errors are positive
negerr = where((*SignalLog.err)[gp] le 0, cnt, COMPLEM=poserr)
if cnt eq n_elements((*SignalLog.err)[gp]) then \$
message, 'The noise is negative or null!!!'
if cnt gt 0 then \$
(*SignalLog.err)[gp[negerr]] = min((*SignalLog.err)[gp[poserr]])

; Check if some pixels have exagerated weight
weight = 1D / ((*SignalLog.err)[gp])^2
large_weight = where(weight gt 100*mean(weight), cnt)
if cnt gt 0 then message, /INFO, \$
'Some pixels have more than 100 times the average weight ... '+ \$
'it may be an error (up to ' + strtrim(max(weight)/mean(weight),2)+')'
endif

;%%%%%%%%%%%%%%%%%%%%  PART III %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%%%%%%%%%%%%%% INITIALIZATION OF THE MODEL TO FIT %%%%%%%%%%%%%%%%%%%%%%%%%%%

lamrange = uly_spect_get(SignalLog, /WAVERANGE, STATUS=status)
if status ne 0 then message, 'Could not obtain WAVERANGE for input spectrum'

velscale = SignalLog.step * c

; uly_fit_init read the models and place them in the cmp.eval_model
status = uly_fit_init(cmp, WAVERANGE=lamrange, VELSCALE=velscale, QUIET=quiet)
if status ne 0 then begin
message, 'Fit initialization failed, abort', /CONT
return
endif

if min(cmp.npix) gt (size(*SignalLog.data))[1]+1 then begin
print, 'The number of pix in the model is greater than in Observation'
print, '* THIS IS PROBABLY A BUG IN THE PROGRAM *'
print, cmp.npix, ' vs. ', (size(*SignalLog.data))[1]
endif

; %%%%%%%%%%%%%%%%%%%  PART IV %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; %%%%%%%%%%%%%%% PREPARE & CALL THE FITTING PROCEDURE %%%%%%%%%%%%%%%%%%
;

; The signal spectrum may be longer than the model, lets cut it
;    (have to cut the obs narrower in case the 2 grids are shifted)
model_range = exp(cmp[0].start + [0.5,cmp[0].npix-1.5] * cmp[0].step)
SignalLog = uly_spect_extract(SignalLog, WAVERANGE=model_range, /OVER)

; initial guess and limits for the LOSVD
;   Note: we have found cases where the convergence in sigma is slow,
;   and where it does not converge. The default of 1 px may not be the
;   best choice. It is OK if sigma < 2 px, but not if sig ~ 10 pix ...
;   Using a large value for the guess seems to converge better.
;   We do not implement it because it is not enough tested
if n_elements(sigma_guess) eq 0  then sigma_guess = SignalLog.step * c  ; 1pix
if n_elements(cz_guess) eq 0  then cz_guess = 0. ; a priori always 0
kguess = [cz_guess, sigma_guess]
if kmoment gt 2 then kguess = [kguess, replicate(0,kmoment-2)]
if n_elements(sigma_limits) eq 2 then klim = [kguess[0] + [-2d3,2d3], sigma_limits]

if not keyword_set(quiet) then begin
velscale = SignalLog.step * c
lamrange = exp(SignalLog.start+[0d,(n_elements(*SignalLog.data))*SignalLog.step])
print,'--------------------------------------------------------------------'
print,'PARAMETERS PASSED TO ULY_FIT'
print,'--------------------------------------------------------------------'
print,'Wavelength range used                 :', lamRange[0], lamRange[1], \$
' [Angstrom]'
print,'Sampling in log wavelength            :', velscale, ' [km/s]'
print,'Number of independent pixels in signal:', \$
ceil(n_elements(*SignalLog.data) / SignalLog.dof_factor)
print,'Number of pixels fitted               :', n_elements(*SignalLog.data)
print,'DOF factor                            :', SignalLog.dof_factor
print,'--------------------------------------------------------------------'
endif

; prepare header for the output structure/file

nloop = 0    ; by default do not loop
if keyword_set(nsimul) then begin
nloop = nsimul - 1
galaxy_log0 = *signalLog.data ; save the input data
endif
if nloop gt 0 then dummy = temporary(file_fits)

for nn=0,nloop do begin

if keyword_set(nsimul) then begin
;       generate a gaussian vector noise for the independent pixels
noisi = randomn(seed,n_elements(galaxy_log0)/signalLog.dof_factor)
;       rebin it as the analysed vector (non-independent pixels)
noise = interpol(noisi,n_elements(galaxy_log0))
noise /= sqrt(mean(noise^2)*signalLog.dof_factor)

noise = randomn(seed,n_elements(galaxy_log0))

*signalLog.data = galaxy_log0 + *signalLog.err * noise
endif

time0 = systime(1)

solution = uly_fit(signalLog,                                     \$
CMP=cmp,                                       \$
KMOMENT=kmoment, KGUESS=kguess,                \$
KFIX=kfix, KLIM=klim, KPEN=kpen,               \$
QUIET=quiet,                                   \$
POLPEN=polpen,                                 \$
CLEAN=clean,                                   \$
MODECVG=modecvg,                               \$
ALL_SOLUTIONS=all_solutions)

if size(solution, /TYPE) ne 8 then begin
if not keyword_set(quiet) then message, /CONT, 'The fit failed, return'
return
endif

; %%%%%%%%%%%%%%%%%%%%%% PART  V %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; %%%%%%%%%%%%%%% PRINT, PLOT AND SAVE THE RESULTS  %%%%%%%%%%%%%%%%%%%%
;

;   make a fits file, which contains the results of the fit
if n_elements(file_fits) eq 1 then begin
uly_solut_swrite, solution, file_fits
endif

;   write the results in the output file
if n_elements(file_res) eq 1 then begin
else uly_solut_twrite, solution, file_res
endif

;   plot final data vs model comparison if required.
if keyword_set(plot) then begin
if n_elements(nsimul) gt 0 then begin
if (nn+1)/100 eq (nn+1)/100. or nn eq nloop then begin
n_para = n_elements(sol[0].losvd)
for n1=0,n_elements(sol[0].cmp)-1 do \$
if ptr_valid(sol[0].cmp[n1].para) then \$
n_para += 1 + n_elements(*sol[0].cmp[n1].para)
n_plot =  (n_para*(n_para-1))/2
wxsize=!D.X_SIZE
wysize=!D.Y_SIZE
plsize = sqrt(wxsize*wysize/n_plot)
nx = floor(wxsize/plsize) ; number of plots in X
ny = floor(wysize/plsize) ; number of plots in Y
miss = n_plot - nx*ny
if miss gt 0 then begin
if ny le nx and miss le ny then nx++ \$
else begin
ny++
if miss gt nx then nx++
endelse
endif
p_save = !P
x_save = !X
y_save = !Y
!P.MULTI = [0, nx, ny]
for n1=1,n_para-1 do for n2=n1+1,n_para do begin
if not keyword_set(quiet) then print, 'plot',n1,n2
uly_solut_tplot, sol, XAXIS=n1, yaxis=n2, QUIET=quiet
endfor
heap_free, sol
!P = p_save
!X = x_save
!Y = y_save
endif
endif else uly_solut_splot, solution
endif

if not keyword_set(quiet) then print, 'time=', systime(1)-time0

if not keyword_set(quiet) then uly_solut_tprint, solution

if not arg_present(solution) then begin
ptr_free, solution.hdr, solution.data, solution.err
heap_free, solution.cmp
endif
endfor

uly_spect_free, SignalLog

END

;--- end ----------------------------------------------------------------------
```
Part of ULySS package