# File rvsao/Emvel/emlfit.x # April 22, 1994 # By Doug Mink, Harvard-Smithsonian Center for Astrophysics # After Bill Wyatt and John Tonry # Copyright(c) 1994 Smithsonian Astrophysical Observatory # You may do anything you like with this file except remove this copyright. # The Smithsonian Astrophysical Observatory makes no representations about # the suitability of this software for any purpose. It is provided "as is" # without express or implied warranty. # EMISSION LINE FITTER include "emv.h" # Note on returned EMPARAMS array - # The elements indexed as [1,2,3],1 are the continuum polynomial # coefficients--there is no midpoint subtracted from the independent var. # The elements indexed as [4,5,6],1 are the fit Gaussian center # (in pixels), sigma (in pixels), and height (in counts). # The elements indexed as [7,8,9],1 are the Equivalent width, Chi**2, # and velocity in km/sec (which is not actually computed in this routine). # The elements indexed as [1-9],2 are the errors in the above, except # element[8],2, which is the number of degrees of freedom. procedure emlfit (fdata, fcont, sky, npix, debug) real fdata[ARB] # Spectrum array with continuum removed real fcont[ARB] # Continuum array real sky[ARB] # Sky spectrum (if non-NULL) int npix # number of pixels bool debug # print diagnostic messages if true int i, j, k, m int lp, rp int ierr # MINI verbosity flag and error return double wlfit # Wavelength to fit on either side of peak(s) int iscombo() int pixleft, pixright, ngauss, nparams, npoly double gcenter[3] # centers for multiple gaussians double gheight[3] # heights for multiple gaussians double gwidth[3] # widths for multiple gaussians double fitvec[20] # working space for 3 Gaussians + 5 cont. params double covar[20] double cont, area, skylevel double wv_cen, wv_cen_plusone double px_cen, px_cen_plusone double dadp, pi int jj, ic, icenter, n int ngfit # Number of Gaussians to fit bool igauss[3] # Flags for found lines in combination real hheight # Half-height of gaussian, cfd, open() double wlj,eqw double wcs_p2w(), wcs_w2p() include "emv.com" begin ierr = 0 pi = 3.14159265358979323846d0 npoly = nlcont # for each line or set of lines, fit the Gaussians i = 1 m = 0 while (i <= nfound) { # See if the line is part of a combination ngauss = iscombo (wlrest, i, nfound, ic, wlfit, ngfit, igauss) # Set number of terms to fit nparams = (ngfit * 3) + npoly if (debug) { call printf ("EMLFIT: line %d %8.3fa: ngauss= %d/%d np= %d\n") call pargi (i) call pargd (wlrest[i]) call pargi (ngauss) call pargi (ngfit) call pargi (nparams) } # Set pixel limits to the fit pixleft = idnint (wcs_w2p (wlobs[i] - wlfit)) if (pixleft < 1) pixleft = 1 pixright = idnint (wcs_w2p (wlobs[i+ngauss-1] + wlfit)) if (pixright > npix) pixright = npix if (debug) { call printf ("EMLFIT: fitting pixels %d to %d\n") call pargi (pixleft) call pargi (pixright) } # Set centers, heights, and widths of found gaussians to fit n = i - 1 do j = 1, ngfit { if (igauss[j]) { n = n + 1 gcenter[j] = pxobs[n] gheight[j] = fdata[int (gcenter[j])] hheight = 0.5 * gheight[j] # Find approx. FWHM of line lp = idnint (gcenter[j]) do k = lp - 1, pixleft { lp = lp - 1 if (fdata[lp] < hheight) break } rp = idnint (gcenter[j]) do k = rp + 1, pixright { rp = rp + 1 if (fdata[rp] < hheight) break } gwidth[j] = double (rp - lp - 1) if (gwidth[j] < 1.d0) gwidth[j] = 1.d0 if (debug) { call printf ("EMLFIT: fit %d %d/%d %8g %8g %8g\n") call pargi (i) call pargi (j) call pargi (ngfit) call pargd (gcenter[j]) call pargd (gheight[j]) call pargd (gwidth[j]) } } } # Set centers, heights, and widths of merged gaussians to fit do j = 1, ngfit { if (!igauss[j]) { do k = 1, ngfit { if (igauss[k]) { gheight[j] = gheight[k] * eht[j,ic] / eht[k,ic] gwidth[j] = gwidth[k] wlj = elines[j,ic] * wlobs[i] / wlrest[i] gcenter[j] = wcs_w2p (wlj) if (debug) { call printf ("EMLFIT: add %d using %d %8g %8g %8g\n") call pargi (j) call pargi (k) call pargd (gcenter[j]) call pargd (gheight[j]) call pargd (gwidth[j]) } } } } } # Fit multiple gaussians ierr = 0 call minifit (fdata, npix, ngfit, gcenter, gheight, gwidth, nparams, pixleft, pixright, debug, fitvec, covar, ierr) # Transfer the continuum poly and each gaussian do j = 1, ngfit { if (!igauss[j]) { next } m = m + 1 do k = 1, 10 { emparams[k,1,m] = 0. emparams[k,2,m] = 0. } # Error? 0=OK, 1=iter. exceeded, 2=Singular if (ierr < 2) { jj = npoly + ((j - 1) * 3) # Continuum polynomial and formal errors if (npoly > 0) { do k = 1, npoly { emparams[k,1,m] = fitvec[k] emparams[k,2,m] = covar[k] } } # Gaussian and formal errors do k = 1, 3 { emparams[k+3,1,m] = fitvec[jj + k] emparams[k+3,2,m] = covar[jj + k] } # Equivalent width, in pixels area = sqrt (2.d0 * pi) * fitvec[jj+2] * fitvec[jj+3] icenter = idnint (fitvec[jj+1]) cont = fitvec[1] + fcont[icenter] if (cont > 0) eqw = area / cont else eqw = cont emparams[7,1,m] = eqw # Error in EQW, pixels if (skyspec) { skylevel = 0.0 do k = fitvec[jj+1] - 2, fitvec[jj+1] + 2 { skylevel = skylevel + sky[k] } skylevel = skylevel / 5.0 } else skylevel = cont area = area + 8.0 * (cont + skylevel) * fitvec[jj+2] + (cont + 2.0 * skylevel) / double (pixright-pixleft+1) * ((eqw) ** 2) if (cont > 0) emparams[7,2,m] = sqrt (abs (area)) / cont else emparams[7,2,m] = sqrt (abs (area)) px_cen = emparams[4,1,m] px_cen_plusone = emparams[4,1,m] + 1.d0 wv_cen = wcs_p2w (px_cen) wv_cen_plusone = wcs_p2w (px_cen_plusone) dadp = wv_cen_plusone - wv_cen # convert EQW to Angstroms emparams[7,1,m] = emparams[6,1,m] * dadp # convert DEQW to Angstroms emparams[7,2,m] = emparams[6,2,m] * dadp # Chi^2, degrees of freedom call gauss (fdata,fitvec,debug,emparams[8,1,m]) emparams[8,2,m] = pixright - pixleft + 1 - (ngauss+1) * 3 } } i = i + ngauss } return end # input list MUST BE WAVELENGTH-SORTED! int procedure iscombo (restw, n, limit, icombo, wlfit, ngfit, igauss) double restw[ARB] # Wavelength sorted list int n # Index of restw to test int limit # Size of list int icombo # Which combination is found double wlfit # Wavelength to fit on either side of peak(s) int ngfit # Number of Gaussians to fit bool igauss[3] # Flags for found lines in combination int i, j, k, m, ngauss include "emv.com" begin icombo = 0 ngauss = 0 m = n igauss[1] = FALSE igauss[2] = FALSE igauss[3] = FALSE do i = 1, ncombo { if (ngauss != 0) break # Check for membership in i'th combination do j = 1, numcom[i] { if (ngauss != 0) break # If member of combination, note that fact if (abs (restw[m] - elines[j,i]) < 0.1) { ngauss = 1 igauss[j] = TRUE wlfit = edwl[i] ngfit = numcom[i] icombo = i # See if another line in the combo is in the input list # by testing the next combo element against next restw line. # restw MUST BE SORTED! do k = j+1, ngfit { m = m + 1 if (m > limit) break else if (abs (restw[m] - elines[k,i]) < 0.1) { ngauss = ngauss + 1 igauss[k] = TRUE } else break } } } } # If not a member of a combination, fit line as single gaussian if (ngauss == 0) { ngauss = 1 ngfit = 1 wlfit = wgfit[n] igauss[1] = TRUE } return ngauss end # Dec 4 1991 Use fit continuum array # Dec 6 1991 Fit all lines in combination whether identified or not # Mar 24 1992 Use combination line information from file # Mar 25 1992 Change GWIDTH to 10 # Mar 31 1992 Avoid dividing by cont if it is =0 # Apr 1 1992 Set fit width in wavelength, not pixels # Nov 9 1992 Allow 0-3 polynomial continuum coefficients # Nov 24 1992 Set fit half-width for single lines from emlines.dat file # Jun 8 1993 Move wavelength <-> pixel conversion to subroutines # Apr 22 1994 Drop length aargument from GAUSS call