;############################################################################# ; This is the version 1.3 of the routine modified in April 2015, based on ; original version (v 1.0) from November 2009. ; PROFIT VERTICAL ; ; Copyright (C) 2009+, Rogemar A. Riffel ; E-mail: rogemar_at_ufsm.br ; ; Please visit http://www.ufsm.br/rogemar/software.html for aplications and updates ; ; If you used this routine, please cite the following reference ; (Riffel, R. A., 2010, ApSS, 327, 239). ; You have permission to modify the current routine for non-comercial personal use. ; Please do not delete the above lines! ; ;############################################################################# ; ; NAME: ; PROFIT ; ; PURPOSE: ; Automatically fit emission/absorption line profiles by Gauss-Hermitte ; series or Gaussian curves from Integral Field Unit observations and obtain ; two-dimensional maps for the gaseous kinematics and flux ; distributions (V, sigma, h3, h4, F) ; ; ; CALLING SEQUENCE: ; PROFIT, INP_FITS, OUT_FILE, REST_LAM, LAM_GUESS, SIG_GUESS, SPECTRAL_WINDOW, $ ; /PLOT, /PRINT, /GAUSSIAN, NITER=10 ; ; INPUT PARAMETERS: ; INP_FITS: Input data cube in standard fits format. Its header must contain ; the CRVAL3, CD3_3, and CRPIX3 ; OUT_FILE: Name of the output multi extension fits file containg the solution ; array. Extension [*,*,0] will contain the flux resulted from the fit; ; [*,*,1] => the centroid velocity; ; [*,*,2] => velocity dispersion; ; [*,*,2] => h3 Gauss-Hermite moment; ; [*,*,3] => h4 Gauss-Hermite moment; ; [*,*,4] => flux obtained directly by the integration of the line profile; ; [*,*,5] => Continuum map ; [*,*,6] => StDev of the continuum map ; [*,*,7] => Chi-square map ; REST_LAM: Rest wavelength of the fitted emission/absorption line ; LAM_GUESS: Initial guess for the centroid wavelength of the line profile ; SPECTRAL_WINDOW: Spectral window in Angstrons used to fit the line profile ; Typicaly in the range from 10 - 50 Angstrons ; PLOT: Set this keyword to plot the resulting fit ; PRINT: Set this keyword to print the solution array ; GAUSSIAN: Set this keyword to fit a Gaussian curve, instead of the ; default Gauss-Hermite series. ; MAX_CHI: Maximum Chi-Square in which the resulting parameters may be used as ; initial guesses for the next fit. If the chi-square is higher than ; MAX_CHI than the initial guesses will be used. ; MAX_SHIFT: If the difference between the position of the peak flux and the ; lambda guess is smaller than MAX_SHIFT (in Angstrons) then this ; position is used as initial guess for the central wavelength. ; OUT_FILE_ERROR: Name of the output fits file containing the error in ; velocity [*,*,1], sigma [*,*,2], h3 [*,*,3] and h4 [*,*,4]. ; NITER: Number of mpfitfun iterations to least-square fitting ; LAM_A, LAM_BX, LAM_BY: Used for initial guess for lambda when chi > max_chi. The initial ; guess for lambda is obtained by guess[3]= LAM_A + LAM_BX * (i+1) + LAM_BY * (j+1) ; ; EXAMPLE: ; profitV, "galaxy.fits", "solution.fits", 12570., 12718, 3., 25, /plot, /print ; ; ; ;############################################################################# ;################## DEFINITION OF THE FITTING FUNCTION ####################### ;############################################################################# FUNCTION fungh, x, par ; Defines the parameters to be fitted a = par[0] b = par[1] ampl = par[2] center = par[3] sigma = par[4] hp3 = par[5] hp4 = par[6] ; creates arrays of zero values y = x*0. GH=x*0. R = x*0. w = x*0. alpha = x*0. H3 = x*0. H4 = x*0. pi = !pi ;Computes the Gauss-Hermite expansion function N=size(GH, /n_elements) FOR i=0,(N-1) DO begin R[i] = a + b * x[i] w[i] = ((x[i]-center)/sigma) alpha[i] = 1./sqrt(2*!dpi) * exp(-w[i]^2/2.) H3[i] = 1. / sqrt(6.) * (2.*sqrt(2.) * w[i]^3 - 3*sqrt(2.)*w[i] ) H4[i] = 1./sqrt(24.) * (4.* w[i]^4 - 12.*w[i]^2 + 3) GH[i] = R[i] + ampl*alpha[i]/sigma*(1. + hp3*H3[i] + hp4*H4[i]) endfor return, GH END ;######################################################################## ;###################### MAIN ROUTINE -- PROFIT ########################### ;######################################################################## pro PROFIT, INP_FITS, OUT_FILE, REST_LAM, LAM_GUESS, SIG_GUESS, SPECTRAL_WINDOW, $ PRINT=PRINT, GAUSSIAN=GAUSSIAN, MAX_CHI=MAX_CHI, MAX_SHIFT=MAX_SHIFT, $ OUT_FILE_ERROR=OUT_FILE_ERROR, NITER=NITER, LAM_A=LAM_A, LAM_BX=LAM_BX, $ LAM_BY=LAM_BY, PLOT=PLOT print, "==================================================================" print, "|| PROFIT: a line-PROfile FITting routine ||" print, "|| R. A. Riffel, 2010, ApSS, 327, 239 ||" print, "==================================================================" ;--21678.6-------------- SOME DEFINITIONS -------------------------- c=299792.458d ; light speed pi=!pi N_par = 7. ; Number of free parameters loadct, 12 ; ----------------- PREPARING THE SPECTRA TO FIT --------------------- fits_read, INP_FITS, cube, hdr ; Read the fits data cube norm=median(cube) cube = cube/norm crval3=sxpar(hdr,'CRVAL3') cdelt3=sxpar(hdr,'CD3_3') crpix3=sxpar(hdr,'CRPIX3') dw=cdelt3 n=size(cube[0,0,*],/n_elements) ilam=(findgen(n)+crpix3)*cdelt3+crval3 xmax=size(cube[*,0,0],/n_elements) ymax=size(cube[0,*,0],/n_elements) dimension=float(xmax*ymax) ; Initialize vectors for each measurement vel=MAKE_ARRAY(xmax,ymax) sig=MAKE_ARRAY(xmax,ymax) h3=MAKE_ARRAY(xmax,ymax) h4=MAKE_ARRAY(xmax,ymax) flux=MAKE_ARRAY(xmax,ymax) fluxa=MAKE_ARRAY(xmax,ymax) continuum=MAKE_ARRAY(xmax,ymax) cntstd=MAKE_ARRAY(xmax,ymax) chi=MAKE_ARRAY(xmax,ymax) sol=make_array(xmax,ymax,9) error=make_array(xmax,ymax,5) ; Outras definicoes iguess = [0.3d, 0.d, 0.5d, LAM_GUESS, SIG_GUESS, 0.0d, 0.00d]; guesses = [a, b, ampl, center, sigma, h3, h4] guess=iguess dl = SPECTRAL_WINDOW ; Distancia em angstrons para cada lado do centro da linha a ser ajustada ; Rotina de ajustes RESOLVE_ROUTINE, 'fungh', /IS_FUNCTION parinfo = replicate({value:0.D,fixed:0.D,limited:[0,0], limits:[0.D,0]}, 7) parinfo[5:6].limited[0:1]=1 parinfo[5:6].limits[0]=-0.5 parinfo[5:6].limits[1]=0.5 parinfo[2].limited[0:1]=1 ; limits for flux parinfo[2].limits[0]=0 parinfo[2].limits[1]=1e30 if keyword_set(gaussian) then begin parinfo[5:6].fixed=1 N_par=5. endif iteration=0 for i=0, (xmax-1), 1 do begin for j=0, (ymax-1), 1 do begin norm1=1 parinfo[*].value=guess cnt1=nint((guess[3] - crval3 - dl)/cdelt3) ; define left limit to continuum fit (in pixels) cnt2=nint((guess[3] - crval3 + dl)/cdelt3) ; define right limit to continuum fit (in pixels) lam=ilam[cnt1:cnt2] spec=cube[i,j,cnt1:cnt2] if avg(spec) gt 0.0 then begin norm1=max(spec,msp) ; Finds the max flux value of spec and its corresponding position (msp) spec=spec/norm1 weights=spec^2 if keyword_set(MAX_SHIFT) then begin ; Use the peak of the profile to define the lambda guess if abs(lam[msp]-guess[3]) lt MAX_SHIFT then begin guess[3]=lam[msp] endif endif ; Fitt each spectrum fitGH=mpfitfun('fungh', lam, spec,weights=weights, parinfo=parinfo, yfit=fitel, perror=perror, $ /nan, /quiet, bestnorm=bestnorm, niter=niter) ;Calculation of the V, sigma, h3 and h4 vel[i,j]=(fitGH[3]-REST_LAM)*c/REST_LAM sig[i,j]=abs(fitGH[4])*c/REST_LAM h3[i,j]=fitGH[5] h4[i,j]=fitGH[6] ; Flux from the fitting func=lam*0 conti=lam*0 func = fitGH[2]*(1./sqrt(2*!dpi) * exp(-((lam-fitGH[3])/fitGH[4])^2/2.))/fitGH[4]*(1.+ fitGH[5]*(1./sqrt(6.)*(2.*sqrt(2.)*((lam- fitGH[3])/fitGH[4])^3 - 3*sqrt(2.)*((lam-fitGH[3])/fitGH[4]))) + fitGH[6]*(1./sqrt(24.) * (4.* ((lam-fitGH[3])/fitGH[4])^4 - 12.*((lam-fitGH[3])/fitGH[4])^2 + 3))) conti=fitGH[0]+fitGH[1]*lam flux_int=func*dw flux[i,j]= total(flux_int)*norm*norm1 ; Flux from the integrated spectrum size=size(spec,/n_elements) conti=conti*0+(avg(spec[0:3])+avg(spec[size-4:size-1]))/2. flux_int=(spec-conti)*dw fluxa[i,j]= total(flux_int)*norm*norm1 ; Averaged continum from the fit continuum[i,j]=avg(fitgh[0]+fitgh[1]*lam)*norm*norm1 cntstd[i,j]=stdev(cube[i,j,cnt1-20:cnt1],cube[i,j,cnt2:cnt2+20])*norm ; Calculates the chi-square spec_med=avg(spec) number_of_pixels=size(fitel,/n_elements) Var = 1./(number_of_pixels-1)*total((spec-spec_med)^2) chi[i,j]=total( ((spec-fitel))^2/Var) / (number_of_pixels - N_par) ; Replace each measurement to the solution array sol[i,j,*]=[flux[i,j], vel[i,j], sig[i,j], h3[i,j], h4[i,j], fluxa[i,j], continuum[i,j], cntstd[i,j], chi[i,j] ] ; Estimates the error of each parameter DOF=N_ELEMENTS(lam) - N_ELEMENTS(fitgh) PCERROR = PERROR * SQRT(BESTNORM / DOF) error[i,j,*]=pcerror[2:6] error[i,j,0]=error[i,j,0]*norm*norm1 ;SOME OPTIONAL KEYWORDS if keyword_set(MAX_CHI) then begin ; If use max_chi to attribute new guesses. if (chi[i,j] lt MAX_CHI) then begin guess = fitGH endif else begin guess=iguess if keyword_set(LAM_A) then begin ; If use max_chi to attribute new guesses. guess[3]= LAM_A + LAM_BX * (i+1) + LAM_BY * (j+1) print, "The lambda guess function have been ran" endif endelse endif if keyword_set(plot) then begin erase plot, lam, spec, xrange=[guess[3]-dl,guess[3]+dl],THICK=2, yrange=[min(spec)*.8,max(spec)*1.05], $ xtitle="Wavelength", ytitle="F!D!7k!X!N (Arbitrary Units)", xtype=0, ytype=0 oplot, lam, fitel, linestyle=0, color=200,THICK=2 oplot, lam, spec-fitel+min(spec), linestyle=1, color=40,THICK=3 endif endif if keyword_set(print) then begin print, '-----------------------------------------------------------------------------------------------------' print, ' Flux V Sigma h3 h4 Fi Chi2' print, sol[i,j,0], sol[i,j,1], sol[i,j,2], sol[i,j,3], sol[i,j,4], sol[i,j,5], sol[i,j,6], FORMAT='(10g12.5)' print, '-----------------------------------------------------------------------------------------------------' endif iteration=iteration+1. print, (iteration/dimension*100.)," % DONE!!" endfor endfor fits_write, OUT_FILE, float(sol) if keyword_set(OUT_FILE_ERROR) then begin error[*,*,1:2] = error[*,*,1:2]*c/REST_LAM fits_write, OUT_FILE_ERROR, float(error) endif end