# File rvsao/Emvel/emgauss.x # June 13, 1994 # By Doug Mink, Harvard-Smithsonian Center for Astrophysics # After Bill Wyatt # 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. #--- Gaussian program interfaces to EMMIN procedure minifit (fdata,fcont, npix, ngauss, gcenters, gheights, gwidths, nparams, pixleft, pixright, debug, retvec, covar, ierr) real fdata[ARB] # Spectrum with continuum removed real fcont[ARB] # Continuum for spectrum int npix # Number of pixels in input spectrum double gcenters[ARB] # Pixel centers to fit double gheights[ARB] # Emission line heights to fit double gwidths[ARB] # Emission line widths to fit int ngauss # Number of gaussians to fit (3 MAX) int nparams # Total number of parameters (15 MAX) int pixleft # Fitting limits int pixright # Fitting limits bool debug # If true, print diagnostic messages double retvec[ARB] # RETURNED vector with minimized parameters double covar[ARB] # RETURNED diagonal of covariance matrix int ierr # RETURNED error flag int i, j int npoly # dimension of continuum fit double conv # converts FWHM to sigma double pcoeff[5] # polynomial coefficients bool gmsperr int section[4] int nsect double ht, wd, htmax, mln2 include "emg.com" data ndata/0/, nderiv/0/ begin mln2 = 0.69314718055994530942d0 npoly = nparams - ngauss * 3 conv = sqrt (8.d0 * mln2) # save in common for communication to CHI**2 routine pixl = pixleft pixr = pixright npar = nparams ngau = ngauss call aclrd (retvec,nparams) # compute initial guess vector # Continuum polynomial nsect = 2 section[1] = pixleft section[2] = pixleft + npoly section[3] = pixright - npoly section[4] = pixright if (npoly > 1) call gmspfit (fdata, nsect, section, npoly, pcoeff, gmsperr) else if (npoly > 0) { pcoeff[1] = 0.d0 gmsperr = FALSE } if (debug) { if (gmsperr) call printf ("MINIFIT: %d-term continuum fit error\n") else call printf ("MINIFIT: %d-term continuum fit used\n") call pargi (npoly) do i = 1, npoly { call printf (" c(%d)= %10.6g ") call pargi (i) call pargd (pcoeff[i]) } call printf ("\n") } # if gmspfit doesn't work, guess another way if (gmsperr) { # estimate slope & offset from end pixels if (npoly >= 2) { retvec[npoly-2] = (fdata[pixright] - fdata[pixleft]) / npix retvec[npoly-1] = fdata[pixleft] - retvec[npoly-2] * pixleft } # constant guess from mean of left and right pixels if (npoly == 1) retvec[npoly-1] = (fdata[pixleft] + fdata[pixright]) / 2.0 } # normal - transfer coefficients else if (npoly > 0) { do i = 1, npoly { retvec[npoly-i+1] = pcoeff[i] } } # Gaussians htmax = 0.d0 do i = 1, ngauss { # Evaluate continuum at center of gaussian ht = gheights[i] wd = gwidths[i] / conv if (ht > htmax) htmax = ht # Set initial values for this gaussian J = npoly + 1 + (i-1) * 3 retvec[j] = gcenters[i] # center retvec[j+1] = ht # height retvec[j+2] = wd # sigma = FWHM/conv } if (debug) { call printf ("MINIFIT: %d-term gaussian fit\n") call pargi (nparams) call printf (" a(i) = ") do i = 1, nparams { call printf ("%10.4g ") call pargd (retvec[i]) } call printf ("\n") } # Fit 1-3 gaussians call emmin (fdata,fcont,nparams,retvec,debug,covar,ierr) if (debug) { call printf (" a(i) = ") do i = 1, nparams { call printf ("%10.7g ") call pargd (retvec[i]) } call printf ("\n") call printf ("da(i) = ") do i = 1, nparams { call printf ("%10.7g ") call pargd (covar[i]) } call printf ("\n") } return end # Calculate chi**2, weighted by the data itself # Called by EMMIN # Note the Chi**2 function is : # # chi^2 = SUM { ( [ Xi - f(Xi) ] **2 ) / weight } # # This routine uses the data itself as the weight. procedure gauss (spectrum,continuum,a,debug,vchi) real spectrum[ARB] # Spectrum real continuum[ARB] # Spectrum double a[ARB] # Polynomial continuum and gaussian parameters bool debug # True to print chi^2 double vchi # RETURNED chi^2 int i double valobs, valcal, valdif, valcon, pix, weight double emfit_eval() include "emg.com" begin vchi = 0.d0 do i = pixl, pixr { if (abs (spectrum[i]) > 1.e-36) { pix = double (i) valobs = double (spectrum[i]) valcal = emfit_eval (npar, ngau, a, pix) valdif = valobs - valcal # data with continuum is the weight valcon = double (continuum[i]) # weight = abs (valobs + valcon) weight = 1.d0 vchi = vchi + (valdif * valdif) / weight } # call printf ("GAUSS: pixel %d: obs: %10.7g cont: %10.7g calc: %10.7g diff: %10.7g\n") # call pargi (i) # call pargd (valobs) # call pargd (valcon) # call pargd (valcal) # call pargd (valdif) } # if (debug) { # call printf ("GAUSS: chi2 is %f\n") # call pargd (vchi) # } return end # Generate 1st derivatives of chi**2 based on parameter vector `a', and # return them in vector `da'. # Called by EMMIN # Note: The first partial derivative of the Chi**2 function is: # # d ( chi**2) d f(Xi) # ---------- = -2 * SUM { [ ( Xi - f(Xi) ) ------- ] / weight } # d A d A # # This routine uses the data itself as the weight. procedure dgauss (spectrum, continuum, a, da, d2a) real spectrum[ARB] # Spectrum real continuum[ARB] # Continuum double a[ARB] # Parameters of function double da[ARB] # First derivatives w.r.t. parameters double d2a[16,16] # Second derivatives w.r.t. parameters int i, j, k, jrow int npoly int npix double pix, gauss, wt, weight double valobs, valcal, valdif, valcon double fact1, fact2, dyda[16] double emfit_eval(), gauss_eval() include "emg.com" begin # zero the return vector call aclrd (da,npar) npoly = npar - ngau * 3 npix = pixr - pixl + 1 # Calculate the derivatives and the Chi**2 derivatives over the fit # and the data points do i = pixl, pixr { # skip if zero weight if (abs (spectrum[i]) < 1.e-36) next # Set double precision variables for computation pix = double (i) valobs = double (spectrum[i]) valcal = emfit_eval (npar,ngau,a,pix) valdif = valobs - valcal # Weight is pixel value including continuum valcon = double (continuum[i]) # weight = abs (valcon + valobs) weight = 1.d0 # -2 * (Obs-Fit) / weight fact1 = -valdif / weight # call printf ("DGAUSS: Pixel %5.0f, obs= %12.7g calc= %12.7g fact1= %12.7g\n") # call pargd (pix) # call pargd (valobs) # call pargd (valcal) # call pargd (fact1) # Polynomial part # call printf (" Polynomial:") do j = 1, npoly { dyda[j] = pix ** (npoly-j-1) # call printf (" %12.7g") # call pargd (dyda[j]) } # call printf ("\n") # Up to 3 Gaussian components at this point do j = npoly+1, npar, 3 { gauss = gauss_eval(pix, a[j]) fact2 = (pix - a[j]) / (a[j+2] * a[j+2]) # (X-C)/sig**2 # Center dyda[j] = fact2 * gauss # Height dyda[j+1] = gauss / a[j+1] # Sigma dyda[j+2] = a[j+2] * fact2 * fact2 * gauss # call printf (" Gaussian= %12.7g, c= %9.3f h= %9.4f s= %9.5f\n") # call pargd (gauss) # call pargd (a[j]) # call pargd (a[j+1]) # call pargd (a[j+2]) # call printf (" fact2= %12.7g, dc= %12.7g dh= %12.7g ds= %12.7g\n") # call pargd (fact2) # call pargd (dyda[j]) # call pargd (dyda[j+1]) # call pargd (dyda[j+2]) # call printf (" Gaussian : %12.7g %12.7g %12.7g\n") # call pargd (dyda[j]) # call pargd (dyda[j+1]) # call pargd (dyda[j+2]) } do j = 1, npar { wt = dyda[j] / weight jrow = (j-1) * npar - 1 do k = 1, npar { d2a[j,k] = d2a[j,k] + wt * dyda[k] } da[j] = da[j] + fact1 * dyda[j] } } return end # Evaluate a polynomial + gaussian # order of coefficients in (a), 15 maximum, is # pn, (continuum polynomial), highest-order coefficient # pn-1, (continuum polynomial) # ...etc... # p0, (continuum polynomial) constant coefficient # g1 (center) # g1 (sigma) # g1 (height) # g2 (center) # ...etc... double procedure emfit_eval (npar, ngauss, a, pix) int npar # total number of parameters int ngauss # Number of Gaussians in the fit double a[ARB] # coefficients double pix # indep. var. int i int npoly double val double gauss_eval(), poly_eval() begin npoly = npar - ngauss * 3 val = poly_eval (pix, a, npoly) do i = npoly+1, npar, 3 { val = val + gauss_eval(pix, a[i]) } return val end double procedure gauss_eval (pix, gcoef) double pix # pixel value double gcoef[3] # Gaussian Center, Height, Sigma double arg, result begin arg = (gcoef[1] - pix) / gcoef[3] result = gcoef[2] * exp (-0.5d0 * arg * arg) return result end double procedure poly_eval (pix, pcoeff, n) double pix double pcoeff[ARB] # Stored highest order first int n double val # RETURNED polynomial value int i begin val = 0.d0 if (n <= 0) return do i = 1, n { val = val * pix + pcoeff[i] } return (val) end # GMSPFIT - fit to multiple sections of an input array # Return is true for success, false otherwise procedure gmspfit (fx, nsect, sections, ndim, coeff, gmsperr) real fx[ARB] # input data array int nsect # number of data sections int sections[ARB] # pairs of (begin,end) points to fit int ndim # dimension of fit double coeff[ARB] # polynomial coefficients bool gmsperr # true if fit done, else false int i, j double chisq pointer x, y, sp int count int nfit begin gmsperr = TRUE if (nsect == 0) return nfit = 0 do i = 1, nsect*2, 2 { nfit = nfit + sections[i+1] - sections[i] + 1 } # allocate d.p. space call smark (sp) x = NULL call salloc (x, nfit, TY_DOUBLE) if (x == NULL) { call printf (" GMSPFIT: cannot allocate y\n") call sfree (sp) return } y = NULL call salloc (y, nfit, TY_DOUBLE) if (y == NULL) { call printf (" GMSPFIT: cannot allocate y\n") call sfree (sp) return } count = 0 do i = 1, nsect*2, 2 { do j = sections[i], sections[i+1] { Memd[x+count] = double (j) Memd[y+count] = double (fx[j]) count = count + 1 } } # Call the fitting routine call polfit (Memd[x], Memd[y], nfit, ndim, coeff, chisq) call sfree (sp) if (chisq > 0.0) gmsperr = FALSE return end # Nov 20 1991 Free x and y before malloc'ing them # Dec 4 1991 Eliminate unused npoly1 # Dec 6 1991 Pass initial values as arguments to minifit # instead of computing them; print final values if debug # Nov 19 1992 Allow a polynomial fits from none to 4th order # add one more point to each section for initial polynomial fit # Nov 24 1992 Add errors to debug listing # Apr 13 1994 Change DCS return line after ftnchek # Apr 22 1994 Do not pass spectrum length to GAUSS, DGAUSS, and D2GAUS # Apr 29 1994 Use stack pointer for GMSPFIT work vectors # Apr 29 1994 Add subroutine to free pointers # May 2 1994 Comment out DHH as it always returns 0 # May 20 1994 Set initial background constant as mean of end pixels # May 25 1994 Add parameter limits # May 27 1994 Increase parameter limits # Jun 6 1994 Set error flag to 3 if parameter limts are exceeded # Jun 7 1994 If constant continuum, set initial value to 0 # Jun 13 1994 Use approximation from Press, et al. for 2nd derivative # Jun 13 1994 Add continuum to GAUSS and DGAUSS argument lists