; RESTRICTIONS: These routines are not approved parts of the ULySS ; package. They are not tested and were designed for a ; particular purpose. ;;loop models pro loop_models model_dir = ['mkb1.30', 'mun0.30', 'mun0.80', 'mun1.30', 'mun1.80', $ 'mun2.30', 'mun2.80', 'mun3.30'] for kkk=0,n_elements(model_dir)-1 do begin grid = uly_ssp_read_vazmiles(uly_root+'/models/Vaz_/website/' + $ model_dir[kkk]+'/') uly_ssp_write, grid, uly_root+'/doc/Vaz_SB_mil_'+model_dir[kkk]+'.fits' heap_free, grid endfor end ;; ;------------------------------------------------------------------------------ ;+ ; NAME: AUX_READ_VAZMILES ; ; PURPOSE: Read a grid of Vazdekis models ; ; USAGE: grid = aux_read_vazmiles (model_dir, $ ; WAVERANGE=waverange, VELSCALE=velscale, $ ; AGERANGE=agerange, METALRANGE=metalrange, $ ; 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 is [3000,10000] ; ; 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 ; Myr and metallicities in dex. ; By default agerange=[0, 23000]Myr, metalrange = [-5.0, +3.0]dex. ; For example ; AGERANGE=[100,1500],METALRANGE=[-1.4,+0.2] ; ; 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 a series of FITS files containing the Vazdekis models ; with the Miles or CaII libraries. It is an analogue to ; AUX_READ_PEGASE. See also ULY_SSP_READ. ; ; EXAMPLE: Read and save a grid of Vaz/Miles to a single fits file ; grid = aux_read_vazmiles(uly_root+'/models/Vaz_/Miles') ; uly_ssp_write, grid, uly_root+'/models/Vaz_Miles.fits' ; ; Read and write the CaT library ; grid = aux_read_vazmiles(uly_root+'/models/Vaz_/caii') ; uly_ssp_write, grid, uly_root+'/models/Vaz_CaT.fits' ; ; LIMITATIONS: We hardcoded number of metallicity bins to 6 and number of age ; bins to 100. The program is working fine if the real ; numbers are equal or less than this but otherwise the ; user have to change. The same is true for the wave-, ; age-, metallicity- ranges. ; ;- ; CATEGORY: AUX ;------------------------------------------------------------------------------ function read_ssp_miles, filename, $ WAVERANGE=waverange, VELSCALE=velscale, $ AGERANGE=agerange, METALRANGE=metalrange, $ ORIGSSP=rr0, $ VERBOSE=verbose ; READ_SSP_MILES (internal routine) ; read Vazdekis models one by one compile_opt idl2, hidden on_error, 2 template_grid_miles = uly_ssp_alloc() extdata = intarr(2,/NOZERO) inum = 0 outwave = dblarr(2,/NOZERO) readext = 0 excl_sb = 0 while ((n_elements(extdata) gt 1) and (readext lt 1)) do begin extdata = mrdfits(filename, inum, mh, /SILENT) naxis1 = double(sxpar(mh,'NAXIS1')) crpix1 = double(sxpar(mh,'CRPIX1')) crval1 = double(string(sxpar(mh,'CRVAL1'))) cdelt = double(string(sxpar(mh,'CD1_1'))) if cdelt eq 0 then cdelt = double(string(sxpar(mh,'CDELT1'))) if crpix1 eq 0 then crpix1 = 1 ; assume that the reference pix is the first one if crval1 le 4 then begin sampling = 1 crval1 *= alog(10d) cdelt *= alog(10d) ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is in log10' endif else if crval1 le 9 then begin sampling = 1 ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is in ln' endif else begin sampling = 0 ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is linear in wavelength' endelse outwave[0] = crval1 - (crpix1 - 1d) * cdelt outwave[1] = outwave[0] + cdelt * (naxis1 - 1d) ; rejection test upon wavelengths ; range if sampling eq 1 then waverange1 = alog(waverange) else waverange1 = waverange if max(waverange1) ge outwave[0]+2*cdelt and min(waverange1) le outwave[1]-2*cdelt then begin wl0 = max([min(waverange1),min(outwave)]) wMax = min([max(waverange1),max(outwave)]) endif else begin message, 'intersection between waverange and template too small: STOP' endelse undefine, waverange1 readext = readext+1 inum = inum+1 endwhile rr0 = readfits(filename, hdr, /silent) ; the wl of the first pixel w0 = double(sxpar(hdr,'CRVAL1')) ; step disp = double(string(sxpar(hdr,'CD1_1'))) if disp eq 0 then disp = double(string(sxpar(hdr,'CDELT1'))) if w0 le 4 then begin sampling = 1 w0 *= alog(10d) disp *= alog(10d) ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is in log10' endif else if w0 le 9 then begin sampling = 1 ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is in ln' endif else begin sampling = 0 ; if not keyword_set(quiet) then $ ; print, 'Assume that the sampling is linear in wavelength' endelse nummin = round((double(wl0)-double(w0))/disp) ; first pixel of the template that we'll use if nummin lt 0 then nummin = 0 wlnummin = w0 + disp * nummin nelem = round( ( double(wMax)-double(max([w0,wl0])) )/disp ) rr0 = rr0[nummin:nelem+nummin] MILES_spectrum=uly_spect_alloc(START=wlnummin, STEP=disp, SAMPLING=sampling,$ DATA=rr0) ; log rebinning for the first age in the template: if n_elements(velscale) gt 0 then begin MILES_log_spectrum = uly_spect_logrebin(MILES_spectrum, velscale) if (size(*MILES_log_spectrum.data))[0] ne 1 then $ message, 'output spectrum must be a vector : STOP' (*template_grid_miles.data) = fltarr(n_elements(*MILES_log_spectrum.data)) (*template_grid_miles.data) = *MILES_log_spectrum.data template_grid_miles.start = MILES_log_spectrum.start template_grid_miles.step = MILES_log_spectrum.step template_grid_miles.sampling = MILES_log_spectrum.sampling uly_spect_free, MILES_log_spectrum endif else begin *template_grid_miles.data = *MILES_spectrum.data template_grid_miles.start = MILES_spectrum.start template_grid_miles.step = MILES_spectrum.step template_grid_miles.sampling = MILES_spectrum.sampling endelse uly_spect_free, MILES_spectrum return, template_grid_miles end ;------------------------------------------------------------------------------ function aux_read_vazmiles, model_dir, $ WAVERANGE=waverange, VELSCALE=velscale, $ AGERANGE=agerange, METALRANGE=metalrange, $ QUIET=quiet compile_opt idl2 on_error, 0 ; Not good to set ranges... if not keyword_set(waverange) then waverange = [3000, 10000] ;angstroms if not keyword_set(agerange) then agerange = [0, 23000] ;Myr if not keyword_set(metalrange) then metalrange = [-5.0, +3.0] ;dex template_grid = uly_ssp_alloc() ; Detect Vazdekis models files = file_search(model_dir[0], '*.fits', COUNT=nfiles) ;print, 'LIST Vazdekis',files if (size(files))[0] eq 1 then begin print,'There are', nfiles,' Vazdekis models' nmet = 7. ;assume there are 7 bins of metallicities nage = nfiles/nmet;100 ; assume age bins (due to the inhomogenious grid of CaII we can not simply make nfiles/nmet) met_table=float(replicate(-999,nmet)) ; initialize with impossible met age_table = make_array(nage, /FLOAT) met_t = make_array(nfiles, /NOZERO, /FLOAT) age_t = make_array(nfiles, /NOZERO) ;; find where the age and metallicity, IMF are located in the fits name char='' k_met = strlen(model_dir[0])-1 while char ne 'Z' do begin char = strmid(files[0],k_met,1) k_met += 1 endwhile if char ne 'Z' then message, 'no metallicity within the fits names: STOP' k_age = k_met+5 while char ne 'T' do begin char=strmid(files[0],k_age,1) k_age += 1 endwhile if char ne 'T' then message, 'no age within the fits names: STOP' imf = strmid(files[0], k_met-7, 6) cnt_F=0 cnt_m=0 cnt_a=0 ; read now each file name and check its features for j=0, nfiles-1 do begin rd_met = strmid(files[j],k_met+1,4) rd_met = (strmid(files[j],k_met,1) eq 'm') ? -1.*rd_met:+1.*rd_met ; test if rd_met is within the fixed range: if metalrange[0] le rd_met and rd_met le metalrange[1] then begin seekm=where(met_table eq rd_met) if seekm[0] eq -1 then begin met_table[cnt_m] = rd_met cnt_m=cnt_m+1 endif rd_age=1000.*float(strmid(files[j],k_age,5)) ; test if rd_age is within the fixed ; range if agerange[0] le rd_age and rd_age le agerange[1] then begin seeka=where(age_table eq rd_age) ;print, age_table, rd_age if seeka[0] eq -1 then begin age_table[cnt_a] = rd_age ;fill the ages cnt_a += 1 ; count the ages endif template_grid_miles = read_ssp_miles(files[j], $ WAVERANGE=waverange, $ AGERANGE=agerange, $ METALRANGE=metalrange, $ VELSCALE=velscale, $ ORIGSSP=origssp, $ VERBOSE=quiet $ ) ; do once at the 1st reading: if cnt_F eq 0 then begin data_t = fltarr(n_elements(*template_grid_miles.data),nfiles,/NOZERO) ; grasp the wcs: template_grid.start = template_grid_miles.start template_grid.step = template_grid_miles.step template_grid.sampling = template_grid_miles.sampling endif data_t[*,cnt_F]= (*template_grid_miles.data) met_t[cnt_F] = rd_met age_t[cnt_F] = rd_age ; incrementing then the counter of allowed files cnt_F += 1 endif endif endfor if cnt_m gt 1 then begin templates = fltarr(n_elements(data_t[*,0]),cnt_a,cnt_m) for i = 0, cnt_F-1 do begin ca = 0 while age_table[ca] ne age_t[i] do begin ca=ca+1 endwhile cm = 0 while met_table[cm] ne met_t[i] do begin cm=cm+1 endwhile templates[*,ca,cm] = data_t[*,i] endfor endif if not keyword_set(quiet) then begin print, 'Nb of files=',cnt_F print print,'metallicities (dex)=' print, met_table[0:cnt_m-1] print,'ages (Myr)=' print, age_table[where(age_table ne 0)] endif *template_grid.data = templates *template_grid.o_age = alog(age_table[where(age_table ne 0)]) *template_grid.o_metal = met_table[0:cnt_m-1] ; grasp general data from the 1st readed fits file sz1 = (size(*template_grid.data))[1] uly_ssp_free,template_grid_miles endif ; reorder the grid in metallicity, THIS HAS A MAJOR BUT ; UNDIAGNOSTICIZED EFFECT ON THE ANALYSIS. PhP 2010/07/26 order = sort(*template_grid.o_metal) *template_grid.o_metal = (*template_grid.o_metal)[order] *template_grid.data = (*template_grid.data)[*,*,order] if not keyword_set(quiet) then begin print, 'Ordered metallicities (dex): ', *template_grid.o_metal help, *template_grid.data endif if n_elements(model_dir) eq 2 then begin ; This is a grid with [Mg/Fe] templates = dblarr(sz1, cnt_a, cnt_m) filename2 = file_search(model_dir[1]+'/' + '*.fits', COUNT=nf) if nf ne cnt_F then $ message, 'The enhanced grid must have the same N of metallicities' for j = 0, nf-1 do begin rd_met = strmid(files[j],k_met+1,4) rd_met = (strmid(files[j],k_met,1) eq 'm') ? -1.*rd_met:+1.*rd_met ; test if rd_met is within the fixed range: if metalrange[0] le rd_met and rd_met le metalrange[1] then begin seekm=where(met_table eq rd_met) if seekm[0] eq -1 then begin met_table[cnt_m] = rd_met cnt_m=cnt_m+1 endif rd_age=1000.*float(strmid(files[j],k_age,5)) ; test if rd_age is within the fixed ; range if agerange[0] le rd_age and rd_age le agerange[1] then begin seeka=where(age_table eq rd_age) ;print, age_table, rd_age if seeka[0] eq -1 then begin age_table[cnt_a] = rd_age ;fill the ages cnt_a += 1 ; count the ages endif template_grid_miles = read_ssp_miles(filename2[j], $ WAVERANGE=waverange, $ AGERANGE=agerange, $ METALRANGE=metalrange, $ VELSCALE=velscale, $ ORIGSSP=origssp, $ VERBOSE=quiet $ ) ; do once at the 1st reading: mm = where(rd_met eq *template_grid.o_metal) aa = where(alog(rd_age) eq *template_grid.o_age) templates[*,aa, mm]= (*template_grid_miles.data) endif endif endfor *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 ; no bad pixels in Vazdekis *template_grid.goodpix = indgen((size(*template_grid.data, /DIM))[0]) ; mkhdr, *template_grid.hdr, *template_grid.data ;*template_grid.data *= 2.e+9 *template_grid.data *= 1.e+7 sxaddpar, *template_grid.hdr, 'BUNIT', '1e-07 LSUN/MSUN/0.1nm', 'Unit of the data' sxaddpar, *template_grid.hdr, 'COMMENT', '--------------------------------------------------' sxaddpar, *template_grid.hdr, 'COMMENT', 'Grid distributed with ULySS http://ulyss.univ-lyon1.fr/' sxaddpar, *template_grid.hdr, 'COMMENT', 'Models originate from http://www.iac.es/proyecto/miles/pages/ssp-models.php' sxaddpar, *template_grid.hdr, 'COMMENT', 'Isochrones Padova 2000' sxaddpar, *template_grid.hdr, 'COMMENT', 'IMF '+imf sxaddpar, *template_grid.hdr, 'COMMENT', 'Age min='+strtrim( exp(min(*template_grid.o_age)),2)+' max='+strtrim( exp(max(*template_grid.o_age)),2)+' Myr' sxaddpar, *template_grid.hdr, 'COMMENT', 'Metallicity min='+strtrim( (min(*template_grid.o_metal)),2)+' max='+strtrim( (max(*template_grid.o_metal)),2)+' dex' if n_elements(model_dir) eq 2 then $ sxaddpar, *template_grid.hdr, 'COMMENT', 'Alpha/Fe 0.0 and 0.4 dex ' sxaddpar, *template_grid.hdr, 'ULY_TYPE', 'SSP' return, template_grid end ;-- end -----------------------------------------------------------------------