; The output ELODIE SSPs of PEGASE.HR have wrong WCS ; fix them here before making the grid for ULySS ; pro fix_phr_elo, dirname filenames = file_search(dirname+'*.fits') for i = 0, n_elements(filenames)-1 do begin h = 0 h = headfits(filenames[i]) ; Read primary header sxaddpar, h, 'CD1_1', 0.2d, 'Step' sxaddpar, h, 'CDELT1', 0.2d, 'Step' sxaddpar, h, 'H_WRESOL', 0.573d, '[0.1nm] FWHM wavelength resolution' modfits, filenames[i], 0, h ; Update header only endfor return end ;------------------------------------------------------------------------------ ;+ ; NAME: AUX_READ_PEGASE ; ; RESTRICTIONS: This routine is not approved part of the ULySS ; package. It is not tested and was designed for a ; particular purpose. ; ; PURPOSE: ; Read a grid of Pegase.HR models ; ; USAGE: grid = aux_read_pegase( model_dir, $ ; WAVERANGE=waverange, VELSCALE=velscale, $ ; AGERANGE=agerange, METALRANGE=metalrange, $ ; MGFE=mgfe, QUIET=quiet ) ; ; ARGUMENTS: ; model_dir: string, directory where the SSP .FITS files are located ; To read models with 2 Mg/Fe enhancements, give ; a string array with the two directories. ; KEYWORDS: ; WAVERANGE: Two element array giving the limits of the range to read, ; in Angstrom. For example WAVERANGE=[4000d,4786d]. ; By default the whole data range is read. ; ; VELSCALE: Size of the pixels in km/s. For example VELSCALE=30. ; If VELSCALE is not given, the models are not rebinned, ie. ; they are as stored in the disk file (in linear or ; log scale). ; ; AGERANGE/METALRANGE: Extracted region of the grid. Ages in ; log(Myr) and metallicities in dex. ; By default the whole data range is read. ; For example ; AGERANGE=alog([100,15000]),METALRANGE=[-1.4,+0.7] ; ; MGFE: Mg/Fe ratio. A value to be added in the header ; ; QUIET: Verbosity control ; ; RETURN: A structure containing the readed models with the ; associeted WCS (see also ULY_SSP_ALLOC) ; ; DESCRIPTION: This routine can be used to generate a new grid of models ; by reading the FITS files produced by the program PEGASE_HR. ; ; It is useful when, for example, a new library, or a ; new version becomes available. ; ; The grid read by AUX_READ_PEGASE can be saved to a file ; and afterward used to define a component with ULY_SSP. ; The set of PEGASE_HR grids of models provided with ; ULySS in the directory 'models/' were generated by ; this program. ; ; EXAMPLE: To read the latest version of PHR/Elodie and save it ; grid = aux_read_pegase(uly_root+'/models/PEG_/e32_Salp_SB') ; uly_ssp_write, grid, uly_root+'/models/PHR_E32_test.fits' ; ; TEST: test_models, uly_root+'/models/PHR_E32_test.fits' ; sp = uly_root+'/data/VazMiles_z-0.40t07.94.fits' ; ulyss, sp, uly_ssp(model=uly_root+'/models/PHR_E32_test.fits') ; test_models, uly_root+'/models/PHR_E32_test.fits' ; ; HISTORY: Martin France, 2006, created ; Mina Koleva, 2007/10/18, debugged, upgraded ;- ; CATEGORY: AUX ;------------------------------------------------------------------------------ function read_ssp_phr, filename, $ WAVERANGE=waverange, VELSCALE=velscale, $ AGERANGE=agerange, METALRANGE=metalrange, $ ORIGSSP=rr0, CONVSSP=rr, $ OUTAGES=outages, $ ; log(age) QUIET=quiet, $ STATUS=status, MODEL_DIR=model_dir ; READ_SSP_PHR (internal routine) ; one-by-one processing of each SSP file compile_opt idl2, hidden on_error, 2 status = 0 ; initialize the error status ;; if filename is not set, or if the pointed file does not exist display a ;; message and return with an error if (filename eq "") then begin message, "no file to read (1st arg is empty)" endif if (file_test(filename) eq 0) then begin message, "filename is not valid (1st arg="+filename+")" endif template_grid_phr = uly_ssp_alloc() ename = 'f' inum = 0 outwave = dblarr(2,/NOZERO) extdata = intarr(2,/NOZERO) readext = 0 excl_sb = 0 time0 = systime(1) fits_open, filename, fcb, MESSAGE=errmsg if errmsg ne '' then begin message, "error reading file="+filename endif while n_elements(extdata) gt 1 and readext lt fcb.nextend do begin extdata = mrdfits(filename, inum, hh, /SILENT, STATUS=status) ename = sxpar(hh, 'EXTNAME') if (ename eq 'ETS_CONTINUUM') then begin hdr = hh rr = extdata naxis1 = sxpar(hh,'NAXIS1') crpix1 = sxpar(hh,'CRPIX1') crval1 = double(string(sxpar(hh,'CRVAL1'))) cdelt = double(string(sxpar(hh,'CD1_1'))) if cdelt eq 0 then cdelt = double(sxpar(hh,'CDELT1')) outwave[0] = crval1 - (crpix1 - 1d) * cdelt outwave[1] = outwave[0] + cdelt * (naxis1 - 1d) if outwave[0] and outwave[1] ne 0 then begin sampling = 0 ; rejection test upon wavelengths range if max(waverange) lt outwave[0] or min(waverange) gt outwave[1] then begin message, 'No intersection between waverange and models: STOP' endif ; search the first pixel to extract from the template. ; waverange indicates the range in the observation as 'center of 1st pix' ; so, the minimum wavelength (edge of 1st pix) is: if n_elements(velscale) gt 0 then $ wl0 = exp(alog(waverange[0]) - 0.5 * velscale / 299792.458d) $ else wl0 = waverange[0] - 0.5 * cdelt w0 = min(outwave) ; the start of the first pixel in the templates is: ww0 = w0 - 0.5d * cdelt nummin = round((double(wl0)-ww0)/cdelt) ; 1st pix of template to extr ; nummin = floor((double(wl0)-ww0)/cdelt) ; 1st pix of template to extr ; search the last pixel to extract from the template. if n_elements(velscale) gt 0 then $ wl1 = exp(alog(waverange[1]) + 0.5 * velscale / 299792.458d) $ else wl1 = waverange[1] + 0.5 * cdelt nummax = round((double(wl1)-ww0+cdelt)/cdelt) ; last pix of template ; nummax = ceil((double(wl1)-ww0+cdelt)/cdelt) ; last pix of template nummin = max([0L, nummin]) nummax = min([fix(naxis1-1), nummax]) wlnummin = w0 + cdelt * nummin ; center of the first extracted pixel endif readext = readext+1 endif else if ename eq 'ETS_CONT_WCA' then begin wavelength = extdata[*].bfit good = where(wavelength ge waverange[0] and wavelength le waverange[1], cnt) if cnt eq 0 then message, $ 'No intersection between waverange and models: STOP' nummin = min(good, MAX=nummax) wavelength = wavelength[good] sampling = 2 readext += 1 endif else if(ename eq 'ETS_PARA') then begin imn = where(extdata.AGE ge agerange[0]) imx = where(extdata.AGE le agerange[1]) ages = [range(imn[0],imx[n_elements(imx)-1])] outages = alog(extdata.AGE) ; rejection test upon ages range if max(alog(agerange)) ge min(outages) and min(alog(agerange)) le max(outages) then begin outages = outages[where(outages ge alog(agerange[0]) and outages le alog(agerange[1]))] endif else begin message, 'intersection between agerange and template too small: STOP' endelse ; rejection test upon metallicities ;; compute the metallicity rigorously (mail from Philippe ;; 2012/2/14 11:07 ) ;; 1996A&AS..117..113G: (Z, Y) = (0,0001, 0,23). ;; ;; 1994A&AS..105...39F: (Z, Y) = (0,1, 0,475). ;; ;; 1994A&AS..105...29F: (Z, Y) = (0,004, 0,24) ; (0,008, 0,25). ;; ;; 1994A&AS..104..365F: (Z, Y) = (0,0004, 0,23) ; (0,05, 0,352). ;; ;; 1993A&AS..100..647B: (Z, Y) = (0,02, 0,28). ; ;; supress the old way which doesn't consider Y (He) ; outmetall = alog10(mean(extdata.Zstars,/double)/0.02d) ;; ; converting to dex z = mean(extdata.Zstars, /DOUBLE) zsun = 0.02d ysun = 0.28d case float(z) of 0.0001: y = 0.23d 0.0004: y = 0.23d 0.004: y = 0.24d 0.008: y = 0.25d 0.02: y = 0.28d 0.05: y = 0.352d 0.1: y=0.475d else: message, 'For Z='+ string(z) + ' no value of Y' endcase outmetall = alog10(z/zsun) - alog10((1-y-z)/(1-ysun-zsun)) if metalrange[0] gt outmetall and outmetall gt metalrange[1] then $ excl_sb = 1 readext += 1 endif inum += 1 endwhile free_lun, fcb.unit ; print, 'reading time1=', systime(1)-time0 ; time 1 is small time0 = systime(1) rr = rr[nummin:nummax, ages] * 1e+07 bunit = '1e-07 LSUN/MSUN/0.1nm' if strmid(filename, strlen(model_dir), 5) eq 'Const' then begin print, 'The SF is constant, multiplying the models by 20000./t' for ii = 0, n_elements(outages) -1 do $ rr[*, ages[ii]] *= 20000./exp(outages[ii]) bunit = bunit+'/20/t' endif sxaddpar, hdr, 'BUNIT', bunit, 'Unit of data' ;print, 'reading time1a=', systime(1)-time0 ;time0 = systime(1) if n_elements(velscale) gt 0 or sampling eq 2 then begin ; rebin the 2D array ; The rebinning it taking 60% of the total model reading time if sampling eq 0 then $ PHR_spectrum = uly_spect_alloc(START=wlnummin, STEP=cdelt, DATA=rr, SAM=0) $ else $ PHR_spectrum = uly_spect_alloc(DATA=rr, SAM=sampling, WAVELEN=wavelength) PHR_log_spectrum = uly_spect_logrebin(PHR_spectrum, velscale) *template_grid_phr.data = *PHR_log_spectrum.data template_grid_phr.start = PHR_log_spectrum.start template_grid_phr.step = PHR_log_spectrum.step template_grid_phr.sampling = PHR_log_spectrum.sampling uly_spect_free,PHR_spectrum uly_spect_free,PHR_log_spectrum endif else begin *template_grid_phr.data = rr template_grid_phr.start = wlnummin template_grid_phr.step = cdelt template_grid_phr.sampling = 0 endelse if excl_sb eq 1 then begin ;outages never set to NULL outmetall = 0d template_grid_phr.start = 0 template_grid_phr.step = 0 template_grid_phr.sampling = 0 endif ;print, 'reading time2=', systime(1)-time0 ;time0 = systime(1) *template_grid_phr.o_age = outages *template_grid_phr.o_metal = outmetall *template_grid_phr.hdr = hdr return, template_grid_phr end ;------------------------------------------------------------------------------ function aux_read_pegase, model_dir, $ WAVERANGE=waverange, VELSCALE=velscale, $ AGERANGE=agerange, METALRANGE=metalrange,$ MGFE=mgfe, QUIET=quiet compile_opt idl2 on_error, 0 if not keyword_set(waverange) then waverange = [3000, 9000] ; 0.1nm if not keyword_set(agerange) then agerange = [1, 23000] ; Myr if not keyword_set(metalrange) then metalrange = [-5.0, +3.0] ; dex template_grid = uly_ssp_alloc() ; detect the fits files in the given directory ; print, model_dir[0] + '/' + '*.fits' filename = file_search(model_dir[0] + '/' + '*.fits') if (n_elements(filename) eq 1) then begin if (filename eq "") then message, "model not found ("+model_dir[0]+")" endif nmetall = n_elements(filename) ; nb of files to be readed ; make an array for the output metallicities ; outmetall = make_array(nmetall,/DOUBLE,/NOZERO) template_grid_phr = read_ssp_phr(filename[0], $ WAVERANGE=waverange, $ AGERANGE=agerange, $ METALRANGE=metalrange, $ VELSCALE=velscale, $ ORIGSSP=origssp, $ QUIET=quiet, $ MODEL_DIR=model_dir $ ) ; grasp general data from the 1st readed fits file sz1 = (size(*template_grid_phr.data))[1] outages = *template_grid_phr.o_age nages = n_elements(outages) if template_grid_phr.step ne 0 then begin template_grid.start = template_grid_phr.start template_grid.step = template_grid_phr.step template_grid.sampling = template_grid_phr.sampling *template_grid.hdr = *template_grid_phr.hdr endif outmetall[0] = *template_grid_phr.o_metal if(nmetall gt 1) then begin templates = dblarr(sz1,nages,nmetall) templates[*,*,0] = (*template_grid_phr.data)[*,*] uly_ssp_free, template_grid_phr for i = 1, nmetall-1 do begin template_grid_phr = read_ssp_phr(filename[i], $ WAVERANGE=waverange, $ AGERANGE=agerange, $ METALRANGE=metalrange, $ VELSCALE=velscale, $ ORIGSSP=origssp, $ QUIET=quiet, $ MODEL_DIR=model_dir $ ) if template_grid_phr.step ne 0 then begin template_grid.start = template_grid_phr.start template_grid.step = template_grid_phr.step template_grid.sampling = template_grid_phr.sampling endif outmetall[i] = *template_grid_phr.o_metal templates[*,*,i] = (*template_grid_phr.data)[*,*] uly_ssp_free, template_grid_phr endfor endif *template_grid.data = templates *template_grid.o_age = outages *template_grid.o_metal = outmetall templates = 0 outages = 0 uly_ssp_free, template_grid_phr if n_elements(model_dir) eq 2 then begin ; This is a grid with [Mg/Fe] templates = dblarr(sz1,nages,nmetall) filename2 = file_search(model_dir[1]+'/' + '*.fits') nmetall2 = n_elements(filename) ; nb of files to be readed if nmetall2 ne nmetall then $ message, 'The enhanced grid must have the same N of metallicities' for i = 0, nmetall-1 do begin template_grid_phr = read_ssp_phr(filename2[i], $ WAVERANGE=waverange, $ AGERANGE=agerange, $ METALRANGE=metalrange, $ VELSCALE=velscale, $ ORIGSSP=origssp, $ QUIET=quiet, $ MODEL_DIR=model_dir $ ) outmetall[i] = *template_grid_phr.o_metal outages = *template_grid_phr.o_age templates[*,*,i] = (*template_grid_phr.data)[*,*] uly_ssp_free,template_grid_phr endfor a = where(outages - *template_grid.o_age,n) if n ne 0 then $ message,'ages from both grids must be the same' b = where(outages - *template_grid.o_age,m) if m ne 0 then $ message,'ages from both grids must be the same' *template_grid.o_mgfe = [0.0, 0.4] ; always these 2 enhancement levels form = [size(templates,/DIM), 2] tot = n_elements(templates) *template_grid.data = reform([reform(*template_grid.data, tot), reform(templates,tot)], form) endif template_grid.title = model_dir[0] *template_grid.waverange = waverange ; determine goodpix lmin = [3000.,5905.,6300.] ; some regions of P.HR are bad lmax = [5880.,6260.,9300.] if template_grid.sampling eq 1 then begin lmin = alog(lmin) lmax = alog(lmax) endif wave = template_grid.start + indgen(sz1) * template_grid.step for i = 0, n_elements(lmin) - 1 do begin goodpixels = where((wave ge lmin[i]) and (wave le lmax[i]), cnt) if cnt gt 0 then begin if n_elements(*template_grid.goodpix) gt 0 then $ *template_grid.goodpix = [*template_grid.goodpix, goodpixels] $ else *template_grid.goodpix = goodpixels endif endfor ;; there is something wrong in the first 3 A of Pegase/Elodie ;; models... problably 0s... i do not know.. NEED TO BE FIXED ;*template_grid.goodpix = (*template_grid.goodpix)[15:*] ; compute the 2nd derivative in case we want to pass this grid to ulyss if size(*template_grid.data, /N_DIM) eq 3 then begin ; normal grid *template_grid.d2t = derive_3d(TMPARR=*template_grid.data, GRID=*template_grid.o_age) endif else begin ; grid with Mg/Fe resolution form = size(*template_grid.data,/DIM) tot = n_elements((*template_grid.data)[*,*,*,0]) *template_grid.d2t = [reform(derive_3d(TMPARR=(*template_grid.data)[*,*,*,0], GRID=*template_grid.o_age), tot),$ reform(derive_3d(TMPARR=(*template_grid.data)[*,*,*,1], GRID=*template_grid.o_age),tot)] *template_grid.d2t = reform(*template_grid.d2t, form) endelse ;; fix the header sxaddpar, *template_grid.hdr, 'REFERENC','2004, A&A 425, 881' sxaddpar, *template_grid.hdr, 'REFCODE', '2004A&A...425..881L' sxaddpar, *template_grid.hdr, 'COMMENT', 'Grid distributed with ULySS http://ulyss.univ-lyon1.fr/' if n_elements(mgfe) then sxaddpar, *template_grid.hdr, 'MGFE', mgfe, 'Mg/Fe ratio' remove = [ $ 'COMMENT ---------- DESCRIPTION OF THE FILE STRUCTURE -------------', $ 'COMMENT This file consists in 3 or 4 Header Data Units',$ 'COMMENT The primary HDU contains the spectra, first axis is the wavelength',$ 'COMMENT direction and the second axis is the Age',$ 'COMMENT The HDU named ETS_LINES lists the emission lines and their',$ 'COMMENT equivalent width',$ 'COMMENT The HDU named ETS_PARA gives the characteristics of each spectrum',$ 'COMMENT Age in Myr, total mas in stars...',$ 'COMMENT If the wavelength sampling is not linear, a fourth HDU, named',$ 'COMMENT If present, ETS_CONT_WCA contains the "dispersion relation",',$ 'COMMENT which is the list of wavelengths of each pixels',$ 'COMMENT Stellar metallicity: 0.10000E-03',$ 'COMMENT Initial metallicity: 0.10000E-03' $ ] for i=0, n_elements(remove)-1 do $ *template_grid.hdr = (*template_grid.hdr)[where(strpos(*template_grid.hdr,remove[i]) ne 0)] sxaddpar, *template_grid.hdr, 'COMMENT', 'Metal min='+string(min(*template_grid.o_metal),F='(f5.2)')$ +' max='+string(max(*template_grid.o_metal),F='(f5.2)')+' dex' sxaddpar, *template_grid.hdr, 'COMMENT', 'Age min='+strtrim(min(exp(*template_grid.o_age)),2)+$ ' max='+strtrim(max(exp(*template_grid.o_age)),2)+' Myr' sxaddpar, *template_grid.hdr, 'COMMENT', 'Evolutionary tracks Padova 1994' sxaddpar, *template_grid.hdr, 'ULY_TYPE', 'SSP' return, template_grid end ;-- end -----------------------------------------------------------------------