; lick, galaxy='NGC3873_restframe.txt', fwhm_galaxy=2.76,template='templates_elodie_tmp.txt',FWHM_star=0.48, velocity=254.2 FUNCTION theSameIndices, a, b na = n_elements(a) nb = n_elements(b) smaller = na ge nb ? nb : na check = na ge nb ? b : a if array_equal(check,a) then compare = b else compare = a ab = make_array(smaller, /INT, value=-1) j = 0 for i=0, smaller-1 do begin ind = where(check[i] eq compare,count) if (count ne 0) then begin ab[j]=check[i] & j++ & endif endfor return, ab[where(ab NE -1)] END PRO stat, a na = n_elements(a) print, min(a), max(a), median(a), total(a)/na END PRO lick, galaxy=gal, fwhm_galaxy=fwhm_gal, template=template, fwhm_star=fwhm_star, velocity=velocity, plot=plot, emission=emission ; Read-in deredshifted galaxy spectrum ; 4795 - 5465 Vazdekis @ fwhm=1.8AA -- template.txt ; galaxy = 'NGC3873_restframe.txt' readcol, gal, F='(F,F,F)', $ wave, flux, galivar, /silent galaxy = flux/max(flux) ; Normalize spectrum to avoid numerical issues galivar = galivar/max(flux) gsize = n_elements(galaxy) ; Read-in Worthey's table of indices in vacuum ; http://astro.wsu.edu/worthey/html/index.table.html fmt='I,F,F,F,F,F,F,F,I,A' readcol, 'lickindices.txt', F=fmt, $ num, band1, band2, blue1, blue2, red1, red2, $ FWHM_lick, unit, indexName, comment='##', /silent sizeindex = (SIZE(indexname))[1] airtovac,band1 airtovac,band2 airtovac,blue1 airtovac,blue2 airtovac,red1 airtovac,red2 ; Smooth galaxy/star spec to Lick resolution of 8.2AA @ 5200AA ; assuming for now constant resolution along stellar spectrum ... ; FWHM of sqrt(3^2-1^2) = 2.82 Angstrom,e.g. ; fsmooth = gaussfold(w,f,2.82) ; Example: Vazdekis spectra have a resolution of 1.8A FWHM at 5100A ; this corresponds to sigma ~ 1.8/5100*3e5/2.355 = 45 km/s. ; FWHM_gal = 2.76 ; SDSS has an instrumental resolution FWHM of 2.76A. ; From blue1 to red2 FWHM == FWHM_lick[i] FWHM = 8.2 ; Lick res at 5200AA FWHM_dif = SQRT(FWHM^2 - FWHM_gal^2) sigma = FWHM_dif/2.355/(wave[1:*]-wave) ; Sigma difference in pixels (log-xaxis) sigma = [sigma[0],sigma] ; a gaussian that extends to 4.0 sigma (99.99%) == default for gsmooth galSmooth = gsmooth(wave,galaxy,sigma) ; CONVOL(ivar^2, gaussfilter, /CENTER, /EDGE_TRUNCATE) ;smooth_err = SQRT(smooth_err) ivar = gsmooth(wave,galivar^2,sigma) ;ivar = sqrt(ivar) ; STAR -------------------------------- if keyword_set(template) then begin readcol, template, F='(F,F)', $ swave, sflux, /silent star = sflux/max(sflux) ; Normalize spectrum to avoid numerical issues ssize = n_elements(star) iwave_both = where(wave GE min(swave) AND wave LE max(swave)) wave = wave[iwave_both] galSmooth = galSmooth[iwave_both] ivar = ivar[iwave_both] FWHM_dif = SQRT(FWHM^2 - FWHM_star^2) sigma = FWHM_dif/2.355/(swave[1:*]-swave) ; Sigma difference in pixels (log-xaxis) sigma = [sigma[0],sigma] ; degrade spectrum to Lick resolution ; a gaussian that extends to 4.0 sigma (99.99%) == default for gsmooth starSmooth = gsmooth(swave,star,sigma) ;stat,swave ;swave_new = wave[where(wave GE min(swave) AND wave LE max(swave))] ;star_new=lin(swave,starSmooth,swave_new) if (n_elements(velocity) GE 1) then begin ; Simple gaussian convolution : velocity = 1.8/5100*cspeed/2.355, 1.8=FWHM ; sigma = FWHM/2.355 ; in Angstroms sigma = velocity[0] * 5100 / 3e5 / (swave[1:*]-swave) ; per pixel sigma = [sigma[0],sigma] if (n_elements(velocity) GT 1) then $ star_conv = ghsmooth(swave,starSmooth,sigma=sigma,hermite=velocity[1:*]) else $ star_conv = gsmooth(swave,starSmooth,sigma) ; -------------------------------------- if (keyword_set(plot)) then begin plot,wave,galSmooth,yrange=[0.,1.] oplot,swave,starSmooth,col='FF00FF'x;,xrange=[4800,4900] oplot,swave,star_conv,col='00FFFF'x pause endif ; /plot endif endif lick = dblarr(25) & lickErr = dblarr(25) & nindex = 25 if (keyword_set(template) AND n_elements(velocity) GE 1) then begin lick=[lick,lick,lick] & lickErr =[lickErr,lickErr,lickErr] & nindex = 75 blue1=[blue1,blue1,blue1] & blue2=[blue2,blue2,blue2] red1=[red1,red1,red1] & red2 = [red2,red2,red2] band1=[band1,band1,band1] & band2 = [band2,band2,band2] unit = [unit,unit,unit] & indexName = [indexName,indexName,indexName] corr = dblarr(25) endif for i=0,nindex-1 do begin if (i LE 24) then begin spec = galSmooth & wave = wave endif else if (i GE 24 AND i LE 49) then begin spec = starSmooth & wave = swave ivar = replicate(1.,n_elements(starSmooth)) endif else begin spec = star_conv & wave = swave ivar = replicate(1.,n_elements(starSmooth)) endelse ; Check spectral range and calculate indices using int_tabulated if possible blue = where(wave GE blue1[i] AND wave LE blue2[i],bcount) red = where(wave GE red1[i] AND wave LE red2[i],rcount) line = where(wave GE band1[i] AND wave LE band2[i],lcount) test1 = where(wave LE blue1[i],testCount1) test2 = where(wave GE red2[i],testCount2) ; If index is covered ; -------------------------------------------------------------- if (testCount1 NE 0 AND testCount2 NE 0) then begin blueMid = (blue1[i]+blue2[i])/2 redMid = (red1[i]+red2[i])/2 centerMid = (band1[i]+band2[i])/2 ; Find pseudocontinuum inside each blue and red interval ; A pseudocontinuum is defined by 'drawing' a line between the ; wavelength midpoints and average flux values of the blue and red ; sides. ; Include fraction of pixels at ends fb1 = lin([wave[blue[0]-1],wave[blue[0]]],[spec[blue[0]-1],spec[blue[0]]],blue1[i]) fb1err = lin([wave[blue[0]-1],wave[blue[0]]],[ivar[blue[0]-1],ivar[blue[0]]],blue1[i]) fb2 = lin([wave[blue[bcount-1]],wave[blue[bcount-1]+1]],[spec[blue[bcount-1]],spec[blue[bcount-1]+1]],blue2[i]) fb2err = lin([wave[blue[bcount-1]],wave[blue[bcount-1]+1]],[ivar[blue[bcount-1]],ivar[blue[bcount-1]+1]],blue2[i]) fb = lin([wave[min(where(wave GE blueMid))], wave[min(where(wave GE blueMid))+1]],[spec[min(where(wave GE blueMid))],spec[min(where(wave GE blueMid))+1]],blueMid) fberr = lin([wave[min(where(wave GE blueMid))], wave[min(where(wave GE blueMid))+1]],[ivar[min(where(wave GE blueMid))],ivar[min(where(wave GE blueMid))+1]],blueMid) fr1 = lin([wave[red[0]-1],wave[red[0]]],[spec[red[0]-1],spec[red[0]]],red1[i]) fr1err = lin([wave[red[0]-1],wave[red[0]]],[ivar[red[0]-1],ivar[red[0]]],red1[i]) fr2 = lin([wave[red[rcount-1]],wave[red[rcount-1]+1]],[spec[red[rcount-1]],spec[red[rcount-1]+1]],red2[i]) fr2err = lin([wave[red[rcount-1]],wave[red[rcount-1]+1]],[ivar[red[rcount-1]],ivar[red[rcount-1]+1]],red2[i]) fr = lin([wave[max(where(wave LE redMid))], wave[max(where(wave LE redMid))+1]],[spec[min(where(wave LE redMid))],spec[min(where(wave LE redMid))+1]],redMid) frerr = lin([wave[max(where(wave LE redMid))], wave[max(where(wave LE redMid))+1]],[ivar[min(where(wave LE redMid))],ivar[min(where(wave LE redMid))+1]],redMid) ; ERRORS as in Cardiel et al. 1998 A&AS 127, 597 disp = [0,wave[1:*]-wave] theta = [wave[line[0]]-band1[i],disp[line],band2[i]-wave[line[lcount-1]]] blueFlux = int_tabulated([blue1[i],wave[blue],blue2[i]],[fb1,spec[blue],fb2])/(blue2[i]-blue1[i]) ;blueFluxErr2 = int_tabulated([blue1[i],wave[blue],blue2[i]],[fb1err,ivar[blue],fb2err]/(blue2[i]-blue1[i])) blueFluxErr2 = 1/(blue2[i]-blue1[i]) * total([fb1err,ivar[blue],fb2err]) ;blueFluxErr2 = total([fb1err,ivar[blue],fb2err]) redFlux = int_tabulated([red1[i],wave[red],red2[i]],[fr1,spec[red],fr2])/(red2[i]-red1[i]) ;redFluxErr2 = int_tabulated([red1[i],wave[red],red2[i]],[fr1err,ivar[red],fr2err]/(red2[i]-red1[i])) redFluxErr2 = 1/(red2[i]-red1[i]) * total([fr1err,ivar[red],fr2err]) ;redFluxErr2 = total([fr1err,ivar[red],fr2err]) if (keyword_set(template) AND n_elements(velocity) GE 1) then begin blueFluxErr2=1. & redFluxErr2=1. & endif sigBlue2 = 1/(blue2[i]-blue1[i])^2 * blueFluxErr2 sigRed2 = 1/(red2[i]-red1[i])^2 * redFluxErr2 iwave = [band1[i],wave[line],band2[i]] ; Continuum flux with fractional pixels contFluxArr = (redMid-iwave)/(redMid-blueMid) * blueFlux + $ (iwave-blueMid)/(redMid-blueMid) * redFlux contFluxErrArr2 = (redMid-iwave)^2/(redMid-blueMid)^2 * sigBlue2 + $ (iwave - blueMid)^2/(redMid-blueMid)^2 * sigRed2 ; Calculate fraction of flux in band1 and band2 points i.e. ; calculate both line (f1,f2) and cont fluxes (c1,c2) and their ; corresponding errors (f1err,f2err,c1err,c2err) in band1 and band2 ; points using lin interpolation f1 = lin([wave[line[0]-1],wave[line[0]]],[spec[line[0]-1],spec[line[0]]],band1[i]) f2 = lin([wave[line[lcount-1]],wave[line[lcount-1]+1]],[spec[line[lcount-1]-1],spec[line[lcount-1]+1]],band2[i]) f1err = lin([wave[line[0]-1],wave[line[0]]],[ivar[line[0]-1],ivar[line[0]]],band1[i]) f2err = lin([wave[line[lcount-1]],wave[line[lcount-1]+1]],[ivar[line[lcount-1]-1],ivar[line[lcount-1]+1]] ,band2[i]) ; Line flux with fractional pixels lineFluxArr = [f1,spec[line],f2] lineFluxErrArr2 = [f1err,ivar[line],f2err] ; Calculate indices using int_tabulated.pro ; lineFlux = int_tabulated(iwave,lineFluxArr/lineFluxErrArr2)/int_tabulated(iwave, lineFluxErrArr2) ; The index itself is an EW, integrating the difference of ; pseudocontinuum and flux over the index passband. ; EW = int[ dlambda * (1-F_lambda/F_cont) ] ; mag = -2.5 * log10(F_lambda/F_cont) = -2.5 * log10(( int[ dlambda * ; ( F_lambda/F_cont) ])/waveRange) ; F_lambda -- flux in the index bandpass; ; F_cont -- flux in the pseudocontinuum (at the center of each index bandpasses ; by lin interpolation in wavelength) ; dlambda -- an index bandpass cov2=0. for k=0, lcount+1 do begin for j=0, lcount+1 do $ if (k NE j) then $ cov2 += lineFluxArr[k]*lineFluxArr[j]/(contFluxArr[k]^2*contFluxArr[j]^2) *$ theta[k]*theta[j] * $ ( (redMid-iwave[k])*(redMid-iwave[j])/(redMid-blueMid)^2*sigBlue2 + $ (iwave[k]-blueMid)*(iwave[j]-blueMid)/(redMid-blueMid)^2*sigRed2) endfor lickErr[i] = sqrt(total((lineFluxErrArr2/contFluxArr^2+lineFluxArr^2*contFluxErrArr2/contFluxArr^4)*theta^2)+cov2) Sc = 1/(band2[i]-band1[i]) * int_tabulated(iwave,lineFluxArr) Cc = (redMid-centerMid)/(redMid-blueMid) * blueFlux + $ (centerMid-blueMid)/(redMid-blueMid) * redFlux sb2 = 1/(blue2[i]-blue1[i])^2 * int_tabulated([blue1[i],wave[blue],blue2[i]],[fb1,spec[blue],fb2]^2) sr2 = 1/(red2[i]-red1[i])^2 * int_tabulated([red1[i],wave[red],red2[i]],[fr1,spec[red],fr2]^2) sig_sb2 = 1/int_tabulated([blue1[i],wave[blue],blue2[i]],[fb1,spec[blue],fb2]^2*[fb1err,ivar[blue],fb2err]) * sb2 sig_sr2 = 1/int_tabulated([red1[i],wave[red],red2[i]],[fr1,spec[red],fr2]^2*[fr1err,ivar[red],fr2err]) * sr2 ; lickErr[i] = Sc/Cc*sqrt(1/int_tabulated(iwave,lineFluxArr^2/lineFluxErrArr2)+$ (redMid-centerMid)^2/(redMid-blueMid)^2 * sig_sb2/Cc^2 + $ ; (centerMid-blueMid)^2/(redMid-blueMid)^2 * sig_sr2/Cc^2) test1 = sqrt(1/int_tabulated(iwave,lineFluxArr^2/lineFluxErrArr2));sqrt(total(lineFluxArr^2*contFluxErrArr2/contFluxArr^4)) test2 = sqrt((centerMid-blueMid)^2/(redMid-blueMid)^2 * sig_sr2/Cc^2);sqrt((redMid-centerMid)^2/(redMid-blueMid)^2 * sig_sb2/Cc^2);sqrt(total(lineFluxErrArr2/contFluxArr^2)) case unit[i] of 0: begin lick[i] = int_tabulated(iwave,1-lineFluxArr/contFluxArr) ;print, indexName[i], lick[i],lickErr[i],0;,Sc/Cc*test1,Sc/Cc*test2, sc,cc end 1: begin lick[i] = -2.5 * alog10(int_tabulated(iwave,lineFluxArr/contFluxArr)/(band2[i]-band1[i])) lickErr[i] = 2.5*alog10(2.73) * 10^(0.4*lick[i]) * 1/(band2[i]-band1[i]) * lickErr[i] ;print, indexName[i],lick[i],lickErr[i],1 end endcase if keyword_set(plot) then begin plot, wave, galaxy, XRANGE=[blue1[i]-20,red2[i]+20], YRANGE=[min(spec[line])-.1,max(spec[line])+.1], title=indexName[i] oplot, [[blue1[i]],[blue2[i]]], [[blueFlux],[blueFlux]],col='FF0000'x, thick=2 oplot, [[blueMid],[blueMid]], [[0],[blueFlux]],col='FF0000'x, line=2;,col=30,line=2 ;col='0000FF'x; oplot, [[red1[i]],[red2[i]]], [[redFlux],[redFlux]],col='0000FF'x, thick=2 oplot, [[redMid],[redMid]], [[0],[redFlux]],col='0000FF'x, line=2;,col=30,line=2 ;col='0000FF'x; oplot,iwave, lineFluxArr, col='00FFFF'x, thick=2 ; xyouts, (band1[i]+band2[i])/2, contFluxArr[0]+0.1, indexName[i],col='00FFFF'x oplot, iwave, contFluxArr, col='00FFFF'x, thick=2 oplot, [[blueMid],[iwave[0]]],[[blueFlux],[contFluxArr[0]]],line=2,thick=2,col='00FFFF'x oplot, [[(iwave)[lcount-1]],[redMid]],[[contFluxArr[lcount-1]],[redFlux]],line=2,thick=2,col='00FFFF'x polyfill, [iwave[0],iwave,iwave[lcount-1]],[contFluxArr[0],lineFluxArr,contFluxArr[lcount-1]],/line_fill,col='00FFFF'x oplot, [[iwave[0]],[iwave[0]]],[[lineFluxArr[0]],[contFluxArr[0]]],thick=2,col='00FFFF'x oplot, [[iwave[lcount-1]],[iwave[lcount-1]]],[[lineFluxArr[lcount-1]],[contFluxArr[lcount-1]]],thick=2,col='00FFFF'x pause endif endif ; -------------------------------------------------------------------------- endfor if (keyword_set(template) AND n_elements(velocity) GE 1) then begin for i=0,24 do begin if (unit[i] EQ 0) then corr[i] = lick[24+i]/lick[49+i] $ ; atomic indices else corr[i] = lick[24+i] - lick[49+i] ; molecular indices print, indexName[i],lick[i],lickErr[i],unit[i], corr[i], round(corr[i]/lick[i] * 100) endfor endif else begin for i=0,24 do $ print, indexName[i],lick[i],lickErr[i],unit[i] endelse ;print,'indices:',lick[0:24] ;print,'corrections:', corr ;print, '-------------------' ;print,corr/lick[0:24] * 100 ; SDSS solution @ ; ui:/home/ana/idlcode/idl/lib/idlspec2d-v5_2_0/pro/science/lickindex/one_indices.pro ; res=8.6/10^wave*0.43429 ;ew[nn]=dw[nn]*(1.0-(ind[nn]/acont[nn])) ;errew[nn]= (dw[nn]*ind[nn]/acont[nn])* $ ; sqrt((errin[nn]/ind[nn])^2+(errcont[nn]/acont[nn])^2) ;magline[nn]=-2.5*alog10(ind[nn]/acont[nn]) ;errmagline[nn]=2.5/2.303*sqrt(errin[nn]^2/ind[nn]^2+errcont[nn]^2/acont[nn]^2) ;SET_PLOT,'PS' ;DEVICE, /COLOR, BITS_PER_PIXEL=8 ;LOADCT, 12 ;angstrom = '!6!sA!r!u!9 %!6!n' ; Lin interpolation between blue and red midpoints ; f = f1 + (f2-f1)/(w2-w1)*(w-w1) ... meaning ; for each w in between (w1,w2) only f(w) > f will be considered for ; integration ; FWHM_tem = 1.8 ; Vazdekis spectra have a resolution FWHM of 1.8A. ; Correct to GH value calculating indices of original and broadened ; stellar spectrum. Correction for ; molecular indices: corr = I_orig - I_broadened ; atomic indices : corr = I_orig/I_broadened ; Read-in best template spectrum ; Calculate indices for this spectrum ; Broaden template with GH(VelDisp,H3,H4) and calculate indices ; Make differences for molecular indices ; Make ratios for atomic indices END