Mike Blanton, Department of Physics, New York University
This web page documents a code written in C and IDL to calculate spectral energy distributions and K-corrections. This version is a significant improvement of an older version still available. Of particular note are that the problem with galaxies at z=0.3 is solved and that the photometric redshift estimates from this version should be more robust.
Standard format IDL documentation exists. Conditions of using the code released here are to send your email address to mb144@nyu.edu, to cite Blanton et al. 2003, and to specify the version tag (eg. v3_2) of the code used. (You can check this in the IDL code using the k_version command).
We document in detail here how to use the SED fits provided to calculate K-corrections and SEDs for particular objects with broad-band magnitude measurements. We do not document in detail how to use the code to optimize the templates; the interested user can peruse the comments in the code and/or contact the authors. Note that although this code has been written with SDSS and 2MASS in mind as its first application, it has been written generally enough to be applied to nearly any set of observations. Although you can trust the results of this version extrapolated to rest-frame wavelengths outside the optical and infrared better than you could the previous versions, extrapolations remain as unwise as they always have been.
The code compiles into a shared object library callable by C, so that people can incorporate the K-correction routines directly into their code. We also provide an interface to the identical library through IDL. Users of Fortran might want to communicate with me about how to accommodate their needs. Or they might want to start using a portable language. It is possible in principle to link the C libraries into code based on SM, TCL/Tk, or Python, and I would be interested in helping interested parties to do this.
One can imagine using the results of this code to calculate the evolution of the luminosity function, the distribution of galaxy colors, as well as to develop galaxy classification algorithms. Or, to estimate photometric redshifts.
First off, if you are downloading the code, please register it. Do so by sending your email address to mb144@nyu.edu. No other information is required. I want this because it will allow me to contact users in case of any improvements or bug fixes.
The best way to obtain the code is to simply download the latest version from this site. The tagged and released versions, from latest to earliest, are:CVS_RSH=ssh export CVS_RSH
cvs -d :pserver:anonymous@spectro.princeton.edu:/usr/local/cvsroot loginWhen the above asks for a password, just hit Enter.
cvs -d :pserver:anonymous@spectro.princeton.edu:/usr/local/cvsroot export -r v3_2 kcorrectwhere you can replace v3_2 with whatever version you want.
It is also possible to CVS "checkout" the latest version of the code. This is not recommended unless you want to become involved in developing kcorrect (that is, if you will want to alter the code and check it back in). Please contact me if you want to do this.
Note that the IDL code depends on the idlutils package.
The code has been tested, and compiles and runs on a Red Hat 7.1 OS. Little experimentation with other operating systems has been performed. Other Linux versions should port very easily, at least. In order to compile the code follow these directions:
KCORRECT_DIR=/usr/local/kcorrect PATH=$KCORRECT_DIR/bin:$PATH LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$KCORRECT_DIR/lib export KCORRECT_DIR export LD_LIBRARY_PATH IDL_PATH=$IDL_PATH:+$KCORRECT_DIR/proor in tcsh:
setenv KCORRECT_DIR /usr/local/kcorrect
setenv PATH $KCORRECT_DIR/bin:$PATH
setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:$KCORRECT_DIR/lib
setenv IDL_PATH ${IDL_PATH}:+$KCORRECT_DIR/pro
These default settings may need to be changed on your system as
follows:
setup -r /usr/local/kcorrect kcorrectThen
Note that the C codes "fit_coeffs.c" and "reconstruct_maggies.c" only depend on the K-correction library; thus, you can use them as templates for incorporating the K-correction code directly into your C code, as long as you include the K-correction library. In fact, you would probably put something like "fit_coeffs.c" into your code but include the k_reconstruct_maggies() call directly inside to get the K-corrections.
We'll describe the basic idea before getting into exactly how to use the code. To each set of five-band galaxy magnitudes I fit an SED which is a nonnegative linear combination of some templates. The templates have been optimized to minimize the residuals between the actual galaxy magnitudes and the galaxy magnitudes reconstructed from the galaxy SED fit. The units of the fit spectra created are in ergs cm^{-2} s^{-1} A^{-1}.
As an example, consider the figure below. I have taken actual spectra from the SDSS and projected them onto the g, r, and i bandpasses (since these are the ones the spectrum overlaps). Then I have run kcorrect and tried to reconstruct the spectrum. What kcorrect does here is to fit the nonnegative sum of three (carefully chosen) templates to the g, r, and i fluxes. The black line is the original spectrum, the red is the reconstructed. As you can see, the fit is reasonable. Not all galaxies are recovered this well. The usual failures are, first, that some galaxies are very red compared to our templates and, second, that the emission lines are not recovered particularly well.
If you want to convert apparent magnitudes in band R to absolute magnitudes in band Q, you need to calculate the K-correction, which is defined by the equation:
m_R = M_Q + DM(z) + K_QR(z)where m_R is the apparent magnitude, M_Q is the absolute magnitude, DM(z) is the distance modulus, accounting for the angular diameter distance and cosmological surface-brightness dimming, and K_QR(z) is the K-correction. By absolute magnitude we mean: the apparent magnitude in band Q the object would have if it were observed at rest, 10 pc away, using an aperture that contains its total flux. See Hogg et al. (2002) for a full description of what the K-correction is. Essentially, the K-correction is fully determined by the galaxy SED and the description of the bandpasses, so that with the SED fits above we have everything we need to calculate K-corrections.
I designed the code to calculate K-corrections between some set of observed AB bands (by default, the SDSS ugriz bands) and those same bands shifted by an arbitrary redshift. For example, for the SDSS spectroscopic survey, most galaxies are observed near redshift 0.1; it makes sense, therefore, to K-correct to bandpasses shifted by 0.1 (ie. by a factor 1.1). This procedure minimizes errors in the estimated K-corrections, because the galaxies at z=0.1 will then have a K-correction which independent of their SEDs (and therefore their colors) and equal to -2.5 log_10 (1+0.1). In general, for a bandpass shifted by z_0, the K-correction for a galaxy at z_0 to the observed bandpass shifted by z_0 is -2.5 log_10 (1+z_0), independent of galaxy SED. (Note that the K-correction is not zero in this case; very early versions of code, v1_10 and previous, incorrectly did not apply this overall offset --- see the note at bottom).
The code is generally based on AB maggies. Two notes:
But you may still be wondering what I mean by an "AB" magnitude. The AB system is designed such that a flat spectrum object with f_nu = 3631 Jy = 3.631e-20 ergs cm^{-2} s^{-1} Hz^{-1} should have every magnitude equal to zero. The beauty of the AB system is that the uniform definition makes it convenient to synthesize AB magnitudes from theory or models. The tragic flaw is that the quality of the AB-ness of a system is very dependent on precise spectrophotometry of standards and the carefulness of the calibrators, since no objects have a flat spectrum. There is a tension between these two needs --- similar to other tensions throughout astronomy between making precise measurements and making interpretable ones.
kcorrect, maggies, maggies_ivar, redshift, kcorrect [ , $ band_shift=, /magnitude, /stddev, lfile=, $ vfile=, vpath=, filterlist=, filterpath=, rmatrix=, $ zvals=, lambda=, vmatrix=, coeffs=, /verbose, /sdssfix ]By default, you hand the code AB maggies in SDSS ugriz, the inverse variance of the maggies (1/stddev^2), and the redshift, and it passes back K-corrections (in magnitudes) to redshift zero for those bandpasses, based on the nonnegative fit of three carefully chosen templates. The flags mean the following:
Please note that the calling procedure has changed since v1_16, in the sense that inputting maggies is now the default. I apologize that some people may have to change their code, but this system is better than magnitudes. Sorry, but then, you get what you pay for.
To perform K-corrections for other sets of bandpasses (say, the Bessell bandpasses or 2MASS bandpasses), look for them in $KCORRECT_DIR/data/filters and insert them in the filterlist= flag.
As an example of using kcorrect, one might do the following:
kcorrect, [1., 4.78, 10.96, 14.45, 19.05], $
[1100., 28., 7.7, 4.4, 2.5], $
0.03, kcorrect, band_shift=0.1, chi2=chi2
which would take the set of maggies in the first argument, observed
for a galaxy at redshift 0.03, and return the K-corrections to the
ugriz magnitudes shifted by 0.1. This choice of maggies and
uncertainties should return a good chi^2. Note that it is a
nonnegative fit, so linear goodness of fit measures don't exactly
apply.
Please note that kcorrect uses a lot of overhead (2-3 seconds worth), so if you can, use it on several thousand galaxies simultaneously at least. If you use it on each object individually (as in the example given above), you will be in trouble. If you truly cannot avoid it, then specify the extra arguments rmatrix and zvals. For example, if you run the following:
kcorrect, [1., 4.78, 10.96, 14.45, 19.05], $
[1100., 28., 7.7, 4.4, 2.5], $
0.03, kcorrect, band_shift=0.1, chi2=chi2, $
rmatrix=rmatrix, zvals=zvals
kcorrect, [1., 4.73, 11.26, 14.25, 18.85], $
[1100., 28., 7.7, 4.4, 2.5], $
0.03, kcorrect, band_shift=0.1, chi2=chi2, $
rmatrix=rmatrix, zvals=zvals
You will find that the second call takes far less time than the first,
because the temporary data stored in rmatrix and zvals will not have
to be regenerated.
If you want to look at the fit spectrum, then specify vmatrix, lambda, and coeffs. The commands
kcorrect, [1., 4.78, 10.96, 14.45, 19.05], $
[1100., 28., 7.7, 4.4, 2.5], $
0.03, kcorrect, band_shift=0.1, chi2=chi2, $
vmatrix=vmatrix, lambda=lambda, coeffs=coeffs
plot, lambda, vmatrix#coeffs, xra=[2000., 12000.]
will plot the best fit spectrum.
If you want to project that spectrum onto a new set of bandpasses, use k_reconstruct_maggies. For example, to project back onto the SDSS bandpasses, the command
k_reconstruct_maggies, coeffs, 0.03, maggies, vmatrix=vmatrix, $
lambda=lambda, filterlist=['sdss_u0.par', 'sdss_g0.par', $
'sdss_r0.par', 'sdss_i0.par', 'sdss_z0.par']
will pass back the reconstruction in the array maggies. You can
project onto any bandpasses you want. For example,
k_reconstruct_maggies, coeffs, 0.03, maggies, vmatrix=vmatrix, $
lambda=lambda, filterlist=['bessell_B.par', 'bessell_V.par']
vega2ab=k_vega2ab(filterlist=['bessell_B.par', 'bessell_V.par'],/kurucz)
bessellmags=-2.5*alog10(maggies)-vega2ab
BminusV=bessellmags[0]-bessellmags[1]
will yield the AB B-V for the best fit. I won't fully document
k_reconstruct_maggies here; look in the IDL documentation for more
information.
The C code requires a bit more attention to use. If people say they would use it if it were easier, I will improve it. In this section, I will describe the stand-alone C programs which fit for the coefficients and which calculate the reconstructed maggies. In the next section, I will briefly describe how to use the libraries within your own C code.
There are two stand-alone programs, fit_coeffs and reconstruct_maggies. fit_coeffs uses the maggies in each band to calculate the coefficients. reconstruct_maggies uses the list of coefficients (as output by fit_coeffs to reconstruct the magnitudes.
fit_coeffs has a help flag:
apeiron.[kcorrect].7 % fit_coeffs --help Usage: catAs indicated, you would call it like| fit_coeffs [--vfile vfile --lfile lfile --ffile ffile ]
cat maggies.dat | fit_coeffswhere each line of "maggies.dat" has the form:
redshift umaggies gmaggies rmaggies imaggies zmaggies uinvvar ginvvar rinvvar iinvvar zinvvarAll of the quantities in this file are in maggies, as described above. So you have to convert the luptitudes or magnitudes --- and the inverse variances --- into maggies before calling this. (Note that the conversion to the inverse variances from the maggies and the magnitude errors is (maggies*0.4*ln(10)*magerr)^{-2}).
From the set of maggies and errors, the code calculates a set of four coefficients which correspond to its guess about the nature of the galaxy SED. The coefficients are output to standard out in the form:
redshift coeff0 coeff1 coeff2where the redshift is included for reasons which will be clear soon.
Once you have the coefficients, you can calculate the reconstructed maggies for a galaxy with the same SED and same total flux, observed at any redshift through a shifted bandpass, using reconstruct_maggies. The help flag for reconstruct_maggies yields
apeiron.[kcorrect].8 % reconstruct_maggies --help Usage: catSo you can call it like:| reconstruct_maggies [--vfile vfile --lfile lfile --ffile ffile --band-shift band_shift --redshift redshift]
cat coeffs.dat | reconstruct_maggieswhere coeffs.dat is in the same format as the output to fit_coeffs:
redshift coeff0 coeff1 coeff2The optional parameter "band_shift" specifies the bandpasses to reconstruct (band_shift=0., the default, indicates to use the observed frame bandpasses). "redshift" specifies the redshift at which to observe the galaxy (by default, the input redshift is used). The output is in the form:
redshift u_rec g_rec r_rec i_rec z_recin the original maggy units.
To summarize with an example, if you wanted to calculate the K-corrections from the observed bandpasses to rest-frame bandpasses shifted by 0.1, you could call the code as follows:
cat $KCORRECT_DIR/test/sample.dat | fit_coeffs | reconstruct_maggies >! maggies.dat cat $KCORRECT_DIR/test/sample.dat | fit_coeffs | reconstruct_maggies --band-shift 0.1 --redshift 0. >! maggies.z0.datThe K-corrections are then defined by -2.5 log_10(maggies/maggies.z0).
You can use the ffile option to reconstruct_maggies to specify a file which contains a list of the filters to use (by default, it uses the SDSS filters).
As I noted above, I have not worked very hard to make the stand-alone C code or the C library very usable, as generally speaking that would have meant replicating features in IDL and C, and I wasn't sure how many people would use the pure C versions. If there is a need to do this, I might be convinced to -- or if YOU do it, we could include your stuff in the library.
Why would you bother incorporating the K-correction code into your code at all, rather than just calculating the K-corrections once and reading in the results? Well, perhaps you have to calculate 1/Vmax, in which for each object you have to check its K-correction at (nearly) arbitrary redshifts. It is best to calculate the K-corrections on the fly in these cases, so we provide ways of doing this through a shared object libary.
For example, the stand-alone C code uses the "libkcorrect.so" library (in fact, this is *exactly* the same library which is linked into IDL). If you have this library in your LD_LIBRARY_PATH and you include the header file "kcorrect.h" (in $KCORRECT_DIR/include), your own code can call it. The best way to figure out how to use it is to read the examples in fit_coeffs.c and reconstruct_maggies.c, but I'll try to explain the essentials here.
$KCORRECT_DIR/data/templates/vmatrix.default.dat $KCORRECT_DIR/data/templates/lambda.default.datThese files are in a special ASCII format which can be read into a "float*" variable in C using the routine k_load_ascii_table. For example:
k_load_ascii_table(&vmatrix,&ndim,&sizes,vmatrixfile);loads in a set of data from vmatrixfile. The pointer to the data is returned in "vmatrix". The number of dimensions in the data is returned in "ndim". "sizes" is an array of size ndim giving the size of each dimension. The information in each of these files is the following:
$KCORRECT_DIR/data/templates/sdss_filters.datThis list is read by the routine "k_load_filters", which assumes all the filters are in the directory $KCORRECT_DIR/data/filters (or if KCORRECT_DIR is not set, in the current directory). This routine is called as:
k_load_filters(&filter_n,&filter_lambda,&filter_pass,&maxn,&nk,filterfile);"nk" is the number of filters. "filter_n" is the number of points in each filter curve, "filter_lambda" is the wavelength scale of each filter, and "filter_pass" is the response curve of each filter. "maxn" is the maximum of "filter_n". This sets the indexing of "filter_lambda" and "filter_pass", so that to get the ith wavelength of the kth filter you access "filter_lambda[k*maxn+i]".
zvals=(float *) malloc(nz*sizeof(float)); for(i=0;i<where nz=1000, zmin=0., and zmax=1.nz;i++) zvals[i]=zmin+(zmax-zmin)*((float)i+0.5)/(float)nz;
rmatrix=(float *) malloc(nz*nv*nk*sizeof(float));
k_projection_table(rmatrix,nk,nv,vmatrix,lambda,nl,zvals,nz,filter_n,
filter_lambda,filter_pass,band_shift,maxn);
One defines "nv" and "nl" based on the "sizes" returned by the
k_load_ascii_table calls described above.
k_fit_noneg(coeffs, rmatrix, nk, nv, zvals, nz, maggies, maggies_ivar,
redshift, ngalaxy, tolerance, maxiter, niter, chi2,
verbose)
The array coeffs should be of the size (nv*ngalaxy), and the coefficients
are returned in it.
k_reconstruct_maggies(zvals, nz, rmatrix, nk, nv, coeffs, redshift,
reconstruct_maggies, ngalaxy)
where reconstruct_maggies are expressed in maggies.
Clever people will realize that a photometric redshift estimation code can easily be constructed from the parts assembled here. Cleverer people will deduce that in fact this has already been implemented and documented in the distribution here. This photo-z code is significantly better than earlier versions (v1_16 and earlier).
As I noted above, I have not worked very hard to make the stand-alone C code or the C library *very* usable, as generally speaking that would have meant replicating features in IDL and C, and I wasn't sure how many people would use the pure C versions. If there is a need to do this, I might be convinced to -- or if YOU do it, we could include your stuff in the library.
It is of course necessary to have filter responses to do any of this work. It turns out that the SDSS filters do not have the responses they were designed to have, but Mamoru Doi has done the hard work of measuring them. Jim Gunn then took these, ran them through models of the telescope and the atmosphere (using an appropriately scaled model of the Palomar atmospheric model at airmass 1.3). Daniel Eisenstein deserves special mention for the extensive testing he has done on these. We have a repository of these and other curves. For kcorrect we use the SDSS files:
sdss_u0.par sdss_g0.par sdss_r0.par sdss_i0.par sdss_z0.parThe other SDSS files refer to the responses for each camcol.
The files are in a special SDSS format called FTCL described at the SDSS DR2 site. IDL code (in particular, yanny_readone()) exists in idlutils to read in such files.
Using this code on a number of SDSS galaxies, we have estimated the linear relationships between a number of different bandpass systems, all in the AB system. This short document shows all of the relevant results.
ERROR: choldc failed (sum==-1.197184e+06)which screws up the inputs to reconstruct_maggies. If this affects you, change the line in k_choldc.c:
printf("ERROR: choldc failed (sum==%e)\n",sum);
to
fprintf(stderr, "ERROR: choldc failed (sum==%e)\n",sum);