|
|
|
|
|
;+ ; 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>] ; [, ADEGREE=<ad>][, MDEGREE=<md>][, POLPEN=<polpen>] ; [, 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 ; by the different tasks. ; ; 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 ; made adding a Gaussian noise equivalent to the estimated ; 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. ; ; The following keywords are handled by the reading function ULY_SPECT_READ, ; 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. ; See more information in ULY_FIT. ; ; ADEGREE ; Default -1; Degree of additive polynomial, ; -1: no 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_read(galaxy, VELSCALE=velscale) ; spectrum = uly_spect_losvdconvol(spectrum, 0., 30., 0, 0, /OVER) ; ulyss, spectrum, MODEL=uly_root+'/models/PHR_Elodie31.fits' ; ; 2. More advanced usages ; 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 ; it under the terms of the GNU General Public License as published by ; 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 function uly_cmp_read, model_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 fits_read, file, d, h, /HEADER_ONLY, MESSAGE=mess 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, $ ADEGREE=adegree, $ 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(adegree) ne 0 then command += ', ADEGREE='+strtrim(adegree,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 n_elements(adegree) ne 1 then adegree = -1 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 cmp = uly_cmp_read(model_file) 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 $ print, 'Degree of additive polynomial ',adegree 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, 'Please report it' 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 sxaddpar, header, 'HISTORY', command 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, $ ADEGREE=adegree, MDEGREE=mdegree, $ 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 if nn eq 0 then uly_solut_twrite, solution, file_res, HEAD=header $ 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 sol = uly_solut_tread(file_res) 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 ----------------------------------------------------------------------