; NAME: ; buildgalaxy.pro (Galaxy Builder) ; ; PURPOSE: ; Creates a galaxy SLD (or SED) from a list of fundamental properties. ; ; CALLING SEQUENCE: ; buildgalaxy, freq, Lnu, [properties] ; buildgalspec, freq, Fnu, z=z, [properties] ; ; INPUTS: ; sfr : current star-formation rate ; mass : total stellar mass ; age : age of stellar population (see below). Can be an array to allow multiple populations. ; model : star-formation history, one of the following: ; 'cont' - continuous, constant star-formation from age through today (default) ; 'ssp' - instantaneous burst at age ; 'exp' - exponentially declining from age to today (may not be reliable) ; texp : exponential decline timesale (for 'exp' model) ; tsmooth : smoothing timescale for all SFH evolution (experimental) ; metallicity : stellar metallicity, in Solar units ; av : dust extinction column in A_V ; embedav : dust extinction column in A_V for heavily dust-embedded stars (for SMG-like galaxies) ; embedfrac : fraction of SFR which is heavily dust-embedded ; dustlaw : attenuation law ; can be: MW, SMC, LMC, 080607, 070802, or Calz (default) ; rv : R_V for dust attenuation ; Tdust : dust temperature in Kelvin, for calculating re-radiation (default 40 K) ; gasmetal : gas-phase metallicity, in Solar units, for calculating line emission ; q : ionization parameter, for calculating line emission ; binning : binning factor, relative to BC03 default resolution ; resolution : smoothing factor ; Ltot : total luminosity (input parameter, displaces sfr) ; sethist : specify star-formation history and time-dependent dust fraction ; nonebular : ignore nebular emission (makes fitting behavior simpler, and relevant if emission lines subtracted) ; pinfo : structure containing elements: sfr, tav, age, age2, tau, tsmooth, metal, mass, fembed, ; av, avembed, dustlaw, rv, Tdust, gasmetal, q, L (for buildgalspec ONLY) ; ; OUTPUTS: ; freq : Frequency of the SED. ; Lnu : Luminosity (/Hz) at this frequency. ; Lbol : output total luminosity (to fix as an input, use Ltot) ; Lir : output luminosity from dust reradiation ; hist : structure showing star-formation history and other time-dependent parameters ; lines : list of lines producing nebular emission (note: strong UV/optical lines only) ; lineluminosity : luminosity of these lines ; ; COMMENTS: ; This code generates population-syntheis galaxy SEDs using a small number of numerical input parameters ; to describe the star-formation history and other basic characteristics. ; Useful for generating plots, and is used as part of the SED fitting software. ; ; The specific parameters required as inputs depend on the model. The simplest possible cases are, e.g.: ; IDL> buildgalaxy, freq, Lnu, age=1e9, mass=1e9, av=0.1 ; The above produces a galaxy SLD for mass of 1e9 solar masses with a single-age stellar population of 1 Gyr. ; IDL> buildgalaxy, freq, Lnu, sfr=1.0, mass=1e9, av=0.1 ; The above produces a galaxy SLD for a mass of 1e9 solar masses and a SFR of 1.0 solar masses per year ; since its formation. (i.e., a constant SFH from 1 Gyr ago until today.) ; Somewhat more complicated star-formation histories can be implemented by using all 3 parameters and inputting ; age as a two-element array denoting the maximum ages of the component stellar populations. For example: ; IDL> buildgalaxy, fnu, Lnu, sfr=1.0, age=[1e8, 1e10], mass=1e10, av=0.1 ; This produces a galaxy with an SFR of 1.0 since 100 Myr ago. But between 100 Myr and 10 Gyr, the SFR was ; something different, such that the total stellar mass was 1e10 solar masses. (The code solves the equation ; mass = age_1*SFR + (age_2-age_1)*SFR_2 to calcualte that SFR, SFR_2.) ; Arbitrarily complex star-formation histories can be specified using the keyword sethist, which contains a ; structure with the array of past SFR values at a grid of (nonlinearly-spaced) ages. That structure can be ; generated by another function that interpolates the function properly. function dustlawrv, dustlaw ; Calculates the default R_V value for the given extinction law (for example, 3.1 for Milky Way). n = n_elements(dustlaw) if n eq 1 then rv = 3.0 else rv = fltarr(n_elements(dustlaw))+3.0 for i = 0, n_elements(dustlaw)-1 do begin dlaw = strlowcase(dustlaw[i]) if dlaw eq 'mw' then rv[i] = 3.1 if dlaw eq 'smc' then rv[i] = 2.74 if dlaw eq 'lmc' then rv[i] = 3.2 if dlaw eq 'calz' then rv[i] = 4.0 if dlaw eq '080607' then rv[i] = 4.0 if dlaw eq '070802' then rv[i] = 3.1 endfor return, rv end function diff, arr ; Calculates the spacing between values in an array. d = abs(shift(arr,1)-shift(arr,-1))/2. d[0] = d[1] d[n_elements(arr)-2] = d[n_elements(arr)-1] return, d end function calcnonthermal, freqGHz, sfr=sfr, alpha=alpha, Lsolar=Lsolar ; Calculate the synchrotron continuum emission component from a galaxy. ; Assumed to be proportional only to the star-formation rate (the total galaxy ; luminosity can also be specified based on the FIR-radio correlation). ; Normalized at L-band (1.4 Ghz) according to Murphy et al. 2011. ; Frequency is input in GHz and luminosity density is returned in cgs units. ; can specify Lsolar rather than SFR ; if FIR-radio: 10^11Lsun galaxy has Lnu@1.4 = 10^22.8 W/Hz = 10^29.8 erg/s/Hz if n_elements(alpha) eq 0 then alpha = 0.75 if n_elements(sfr) eq 0 then begin ; use luminosity, via FIR-radio correlation Lnu14 = 10.0^(29.8) * (Lsolar/1e+11) ; Lnu at 1.4 GHz sfr = Lnu14 / (4d+28*0.45) endif else begin ; use current SFR ; ; from Carilli & Yun ;Lnu14 = 4d+28 * sfr*0.45 ; needed to add 0.45 for IMF>SNmass correction ;Lnu14 = 3d+28 * sfr ; based on eq. 13 in Yun&Carilli (note Jy and Mpc), ; not used since they seem to k-correct wrong fth = 1.0 ; Condon: L = 5.3d+21 W/Hz (nu/GHz)^(-alpha) * SFR(>5) ; = 5.3d+28 erg/s/Hz (nu/GHz)^(-alpha) * SFR(>5) ;Lnu1 = 5.3d+28 * fth * sfr * 0.45 ; now working at 1.0 GHz directly following Condon ; ; (for alpha=0.75, would give 4.1, same as Carilli) ; Murphy: Lnu14 = (1./6.35d-29) * SFR endelse ;dLcm = dLMpc * (3.086d+24) ;fcgs = Lnu14 * (1+z)^(1+alpha) / (4*3.1415*dLcm^2) * (freqGHz/1.4)^(-alpha) ;Lnu = Lnu1 * (freqGHz/1.0)^(-alpha) Lnu = Lnu14 * (freqGHz/1.4)^(-alpha) Lnu *= bkpw(freqGHz/5000,-3+alpha,0,sharp=1) ; add a break at high frequencies ;fuJy = fcgs*1d+29 return, Lnu end function calcthermal, freqGHz, T=T, Lsolar=Lsolar, a=a, hotT=hotT, hota=hota, hotfrac=hotfrac, vhotT=vhotT, vhota=vhota, vhotfrac=vhotfrac, beta=beta ; Calculate the thermal continuum due to dust re-radiation. ; The total *dust* luminosity is input as Lsolar. ; Most of the energy (55% by default) is reradiated at (dust) temperature T, which is typically around ; 40 K (default value). ; However, two components of hotter, smaller grains are also included. The temperatures can be specified ; as hotT and vhotT; by default they are 2x and 4x the input temperature. Their fractional luminosities can ; be specified as hotfrac and vhotfrac; by default they are 40% and 5% (i.e., 0.4 and 0.05). ; The characteristic grain size (above which radiation is inefficient) is given in microns by a for the standard ; dust and hota and vhota for the hotter dust. ; The spectral index at wavelengths above the characteristic grain size is given by beta. ; ; A variety of real-world SEDs can be fit by adjusting all of these parameters, although this typically requires excellent ; data. In general, most of them are unknown and default values are employed: typically only T and Lsolar are set. ; Frequency is input in GHz and luminosity density is returned in cgs units. ; Lsun = 3.84e+33 erg/s ; nu_c = 2000. ; 2e12 Hz (Carilli & Yun), equiv. of 2pia=150. ; deprecated way to specify long-wave emission if n_elements(beta) eq 0 then beta = 1.0 if n_elements(a) eq 0 then a = 0.1 ; micron if n_elements(T) eq 0 then T = 40. ; default dust temperature is 40 K if n_elements(hotT) eq 0 then hotT = (T < 50) * 2 if n_elements(hota) eq 0 then hota = a / 2 if n_elements(hotfrac) eq 0 then hotfrac = 0.25 ;0.4 if n_elements(vhotT) eq 0 then vhotT = (T < 50) * 4 if n_elements(vhota) eq 0 then vhota = a / 4 if n_elements(vhotfrac) eq 0 then vhotfrac = 0.025 ;0.05 coolfrac = 1-hotfrac-vhotfrac ; Three-phase dust model replicates the blue side of the SED pretty well. ; There are two stochastically-heated and-or smaller dust components which dominate in the short-wavelength regime. ; note that dust grains are small, so this is NOT a Planck function! ; a = grain size in microns beta = emissivity ; cool dust plancknorm = (freqGHz^3/(exp(4.8d-2*freqGHz/T)-1) / T^4. / 1.2162377d+15) ; total P*dnu is 1 erg/s Qdust = ((2*3.1415*a)/(3e+5/freqGHz))^beta < 1 ;Qdust = 1 - exp(-(freqGHz/nu_c)^beta) cooldustflux = coolfrac * Qdust * plancknorm ; hot dust plancknorm = (freqGHz^3/(exp(4.8d-2*freqGHz/hotT)-1) / hotT^4. / 1.2162377d+15) ; total P*dnu is 1 erg/s Qdust = ((2*3.1415*hota)/(3e+5/freqGHz))^beta < 1 ;Qdust = 1 - exp(-(freqGHz/(2*nu_c))^beta) hotdustflux = hotfrac * Qdust * plancknorm ; very host dust plancknorm = (freqGHz^3/(exp(4.8d-2*freqGHz/vhotT)-1) / vhotT^4. / 1.2162377d+15) ; total P*dnu is 1 erg/s Qdust = ((2*3.1415*vhota)/(3e+5/freqGHz))^beta < 1 ;Qdust = 1 - exp(-(freqGHz/(4*nu_c))^beta) vhotdustflux = vhotfrac * Qdust * plancknorm dustflux = cooldustflux + hotdustflux + vhotdustflux bad = where(finite(dustflux) eq 0, ct) if ct gt 0 then dustflux[bad] = 0 ; not sure how this happened - freq apparently went to infinity... ;dnu = abs(shift(freqGHz,1)-freqGHz)*1d9 ;dnu[0] = dnu[1] dnu = diff(freqGHz)*1d9 ;renorm = total(plancknorm*dnu) / total(dustflux*dnu) ;dustflux = dustflux * renorm Lnucgs = (Lsolar*3.84d+33) * dustflux / total(dustflux*dnu) ;LnuuJy = Lnucgs*1d+29 return, Lnucgs end function calcpah, freqGHz, Lsolar=Lsolar ; Calculate simplified polycyclic aromatic hydrocarbon emission SLD. ; The emission line ratios and widths are specified in a text file and always the same. ; Only the total luminosity varies (Lsolar). ; (In reality, the line ratios do vary and there are non-trivial absorption corrections.) common galbuilder, galbuilderpath common pahblock, pah if n_elements(galbuilderpath) eq 0 then galbuilderpath = getgalbuilderpath() if n_elements(pah) eq 0 then begin readcol, galbuilderpath+'pahlines.txt', pahwav, pahgamma, pahfwhm, pahlum, /silent, comment='#' pah = replicate({wav:0., gamma:0., fwhm:0., lum:0., freq:0.}, n_elements(pahwav)) pah.wav = pahwav pah.gamma = pahgamma pah.lum = pahlum pah.fwhm = pahfwhm pah.freq = 2.998e5 / pah.wav endif npah = n_elements(pah) pahflux = fltarr(n_elements(freqGHz)) d = abs(alog10(freqGHz/25000.)) xfact = fltarr(n_elements(freqGHz))+1 xfact[where(d ge 1.2)] = (1+d[where(d ge 1.2)]-1.2)^(-2.) ; PAHs actually contribute significant flux to radio in some circumstances for a Lorentzian, so suppress the outer wings for i = 0, n_elements(pah)-1 do begin pahflux += pah[i].lum * (pah[i].gamma*pah[i].freq/2)^2 / ( (freqGHz-pah[i].freq)^2 + (pah[i].gamma*pah[i].freq/2)^2 ) endfor pahflux *= xfact ;dnu = abs(shift(freqGHz,1)-freqGHz)*1d9 ;dnu[0] = dnu[1] dnu = diff(freqGHz)*1d9 tot = total(pahflux*dnu) Lnucgs = (Lsolar*3.84d+33) * pahflux / tot ; normalize to specified luminosity return, Lnucgs end function redden, wav, flux, av, rvin, dustlaw=dustlaw ; simple dust-reddeining routine if n_elements(rvin) eq 0 then rv=dustlawrv(dustlaw) if n_elements(rvin) gt 0 then if rvin lt 0 then rv=dustlawrv(dustlaw) if n_elements(rv) eq 0 then rv = rvin ; should ultimately use the lext function (set the power-law index to 0 and multiply) ebv = av/rv llaw = strlowcase(dustlaw) if llaw eq 'mw' then fm_unred, wav, flux, -ebv, r_v=rv if llaw eq 'smc' then fm_unred, wav, flux, -ebv, r_v=rv, c1=-4.96, c2=2.26, c3=0.389, c4=0.46 if llaw eq 'lmc' then fm_unred, wav, flux, -ebv, r_v=rv, /avglmc if llaw eq '080607' then fm_unred, wav, flux, -ebv, r_v=rv, c3=1.3, c4=0.3 if llaw eq '070802' then fm_unred, wav, flux, -ebv, r_v=rv, c3=2.68, c4=0.32, gamma=1.09, x0=4.642 if llaw eq 'calz' then calz_unred2, wav, flux, -ebv, r_v=rv return, flux end function calcnebular, wav, sfr0=sfr0, gasmetal=gasmetal, resolution=resolution, q=q, $ lines=lines, linewav=linewav, lineluminosity=lineluminosity ; Calculates SLD for nebular emission lines (UV/optical lines excited by massive stars) ; Same general philosophy as for other calc* functions, but note that this is called from ; within calcstellar and the input is wavelength in Angstroms, *not* frequency. ; (todo: should substitute most of this out with calcnebularluminosity, which is supposedly better) common kewleyblock, kewley if n_elements(kewley) eq 0 then kewley = readkewley() lineflux = fltarr(n_elements(wav)) ;wav = 2.998e9/freqGHz if n_elements(resolution) eq 0 then resolution = 1. x = alog10(gasmetal)+8.69 ; x = [log[O/H]+12], solar is 8.69 (Asplund) if n_elements(q) eq 0 then q = 3d7 ; supposedly typical from Kewley lines = ['Halpha', 'Hbeta', 'Hgamma','Hdelta','Hepsilon','Hzeta', 'Lyalpha',$ '[OII]1','[OII]2', '[OIII]1','[OIII]2', '[NII]1','[NII]2', '[SII]1', '[SII]2'] linewav = [6562.82,4861.33,4340.5,4101.7,3970.1,3889.1, 1217., $ 3726.1,3728.8, 4958.91,5006.84, 6548.06,6583.57, 6716.44,6730.82] linelumh = 1.255d+41*[1., 1./2.86, 1/6.08, 1./(2.86/0.260), 1./(2.86/0.159), 1./(2.86/0.105)] ; http://www.naoj.org/Science/Resources/lines/hi.html , itself from Hummer & Storey (1987; MNRAS 224, 801). ; original reference is probably Kennicutt, Tamblyn, & Congdon (1994) linelumlya = 8.7 * linelumh[0] ; Dijkstra & Westra k = igrid(kewley.oiihalpha, kewley.q, q) linelumoii = 1.265d+41*poly(x,k) * [0.5,0.5] if x lt 8.2 then linelumoii = 1.265d+41*poly(8.2,k) * [0.5,0.5] * 10.^(x-8.2) ; Kewley OII/Halpha goes negative eventually, and the grid doesn't deal with metallicity below 8.2... ; going to need to look into this later...! ; these come from Fig 9 / Table 4 of Kewley 2004 ; New behavior at low-Z is based on roughly matching the OIII/Hbeta plot k = igrid(kewley.logoiiioii, kewley.z, gasmetal) linelumoiii = total(linelumoii) * 10.0^(poly(alog10(q),k)) * [0.33333,1.0] k = igrid(kewley.logniihalpha, kewley.q, q) linelumnii = linelumh[0] * 10.0^(poly(x,k)) * [0.339,1.0] ;k = igrid(kewley.logniioii, kewley.q, q) ;linelumnii = total(linelumoii) * 10.0^(poly(x,k)) * [0.339,1.0] k = igrid(kewley.logniisii, kewley.q, q) linelumsii = linelumnii[1] / 10.0^(poly(x,k)) * [0.6,0.4] linelum = [linelumh, linelumlya, linelumoii, linelumoiii, linelumnii, linelumsii] linelum = linelum*sfr0 / 3.839d+33 ; where does this come from? for i = 0, n_elements(lines)-1 do begin ;near = where(abs(wav - linelum) lt peak = linelum[i]/(resolution*sqrt(2*!pi)) lineflux += gaussian(wav,[peak, linewav[i], resolution, 0]) / 1.0 ; per Ang endfor ; return luminosities now and fluxes get converted later, for some reason lines = ['H_alpha', 'H_beta', 'H_gamma', '[OII]', '[OIII]', '[NII]', 'Ly_alpha'] linewav = [6562.82,4861.33,4340.5, 3727., 5006.84, 6583.57, 1217.] lineluminosity = [linelumh[0], linelumh[1], linelumh[2], total(linelumoii), max(linelumoiii), max(linelumnii), linelumlya]*sfr0 ;print, q, gasmetal, '|', lineluminosity/lineluminosity[0] ; wav is elsewhere, change that too if you change it here!!! return, lineflux end pro calcstellar, freqGHz, Lnu, $ sfr_in=sfr_in, mass_in=mass_in, $ tav_in=tav_in, age_in=age_in, texp_in=texp_in, tsmooth_in=tsmooth_in, metallicity_in=metallicity_in, $ embedfrac_in=embedfrac_in, embedtime_in=embedtime_in, av_in=av_in, embedav_in=embedav_in, $ dustlaw_in=dustlaw_in, rv_in=rv_in, $ gasmetal_in=gasmetal_in, q_in=q_in, $ binning=binning, resolution=resolution, $ ; vissfr=vissfr, $ Lbol=Lbol, Lir=Lir, hist=hist, sethist=sethist, nonebular=nonebular, $ lines=lines, lineluminosity=lineluminosity ; The key subroutine of galaxy builder - calculates the star-formation history and sums the SED out of BC03 components. ; Most inputs are as for buildgalaxy, below. The implementation was designed for a lot of flexibility in specifying multiple ; components but for now only 2 components are allowed, and only age can be vectorized using the keyword parameters. ; Additional inputs that currently cannot be specified at buildgalaxy level: ; embedtime: maximum age for a SED component to be heavily obscured. ; internal native units are Lsol/ang (I think..) but returns freq and Lnu(cgs) ; output defaults are zero in case something fails (invalid input, etc.) and we return early ; freqGHz = 10.0^(3.5+3.5*fltarr(100)/100) ; why was this here? Lnu = [0., 0., 0.] Lir = 0. Lbol = 0. dum = fltarr(221) sfhist = 0*popsynhist(age=1e9, model='ssp') lines = [''] lineluminosity = [0.] Lnu = popsynsed(wav, sfhist, mass=0., bct=bct, $ ; acutally Llambda, doesn't matter for null case binning=binning, resolution=resolution) ; new, hopefully OK ; bizarre gap between 10-30 microns hist = {t:bct, sfr:sfhist, sfrobs:sfhist} freq = 3d+18/wav ; Hz freqGHz = freq / 1e+9 ; already did this if n_elements(mass_in) gt 0 then if mass_in eq 1.0 then begin ;print, 'Mass was set to 1.0 - failed fit???' ; print an error message if the print error messages flag that you will implement soon is set return endif ; first of all, determine how many populations we have. ;npop = max([n_elements(age_in), n_elements(texp_in), n_elements(mass_in)]) npop = 1 if n_elements(age_in) eq 2 then begin if age_in[0] gt 0 and age_in[1] gt 0 then npop = 2 if keyword_set(tav_in) then if tav_in gt 0 then npop = 2 endif ; the second component must have a specified age XXX relaxing this for now ; n_elements(Av_in), n_elements(embedfrac_in), n_elements(embedtime_in), n_elements(embedav_in), n_elements(metallicity_in), ;print, mass_in ;print, age_in ;print, sfr_in ;print, texp_in ;print, npop ; transfer inputs to real variables - this is to be sure we don't actually touch the outputs ; negative values are treated as unset. You cannot specify an array to be partially negative/default. if n_elements(sfr_in) gt 0 then if sfr_in[0] ge 0 then sfr = sfr_in if n_elements(mass_in) gt 0 then if mass_in[0] ge 0 then mass = mass_in if n_elements(tav_in) gt 0 then if tav_in[0] ge 0 then tav = tav_in if n_elements(age_in) gt 0 then if age_in[0] ge 0 then age = age_in if n_elements(texp_in) gt 0 then if texp_in[0] ge 0 then texp = texp_in if n_elements(tsmooth_in) gt 0 then if tsmooth_in[0] ge 0 then tsmooth = tsmooth_in if n_elements(metallicity_in) gt 0 then if metallicity_in[0] gt 0 then metallicity = metallicity_in if n_elements(embedfrac_in) gt 0 then if embedfrac_in[0] ge 0 then embedfrac = embedfrac_in if n_elements(embedtime_in) gt 0 then if embedtime_in[0] gt 0 then embedtime = embedtime_in if n_elements(av_in) gt 0 then if av_in[0] ge 0 then av = av_in if n_elements(embedav_in) gt 0 then if embedav_in[0] ge 0 then embedav = embedav_in if n_elements(dustlaw_in) gt 0 then if strmid(dustlaw_in,0,1) ne '-' then dustlaw = dustlaw_in if n_elements(rv_in) gt 0 then if rv_in[0] gt 0 then rv = rv_in if n_elements(gasmetal_in) gt 0 then if gasmetal_in[0] gt 0 then gasmetal = gasmetal_in if n_elements(q_in) gt 0 then if q_in[0] gt 0 then q = q_in ; Simple defaults if n_elements(texp) eq 0 then texp = !values.f_infinity ;100e+6 if n_elements(tsmooth) eq 0 then tsmooth = 0 if n_elements(metallicity) eq 0 then metallicity = 1.0 if n_elements(embedfrac) eq 0 then embedfrac = 0.5 if n_elements(embedtime) eq 0 then embedtime = 200e+6 if n_elements(av) eq 0 then av=0. ; making this vector is nontrivial if n_elements(embedav) eq 0 then embedav = 50 if n_elements(dustlaw) eq 0 then dustlaw = 'Calz' if n_elements(rv) eq 0 then rv=dustlawrv(dustlaw) if n_elements(gasmetal) eq 0 then gasmetal = metallicity if n_elements(sethist) eq 0 then begin ; bypass all this with a custom history ; Complex defaults : (SFR), mass, age. ; specify sfr, scalar mass, vector age -> controls mass ratio. ;print, sfr, mass, age[0], age[1], sfr*min([age[0],age[1]]) gt mass, npop ; multiple populations with partially specified age if n_elements(sfr) gt 0 and n_elements(mass) eq 1 and npop eq 2 and n_elements(age) eq npop then begin if n_elements(tav) gt 0 then begin m1 = age[0]*sfr ; mass more recent than impulse if m1 gt mass then begin print, 'SFR/age/mass (',sfr,tav,mass,') are incompatible for this model.' return endif m2 = mass - m1 ; mass older than impulse, different from vector mass which is blocked vertically ; this permits age[1] < age[0]. if m2 gt 0 then age[1] = (mass/m2) * (2*tav - age[0]) $ else age[1] = 10*age[0] ; avoid 0/0 if age[1] lt age[0] then begin ;print, mass/m2, 2*tav, age[0], ' ->', age[1] ;print, 'Inverted population ages are not allowed.' return endif ;print, tav, age[0], age[1], m1, m2, mass, sfr endif else begin ;if min(age)*sfr gt mass then print, 'Impossible age*SFR/mass, setting age to mass/SFR' if min(age)*sfr gt mass then return ; age[where(age eq min(age))] = mass/sfr ; effectively make it a 1-population model by ; ignoring input age and setting it to real value. ; this makes chi^2 flat, slows convergence, and ; age parameter will appear meaningless. ; problematic if age was fixed. ;if max(age)*sfr lt mass then print, 'Declining SFH detected, setting age to mass/SFR' ;if max(age)*sfr lt mass then age[where(age eq max(age))] = mass/sfr endelse if finite(texp) then sfrpermass = 1. / (texp * (exp(age/texp) - 1)) $ ; this is a vector else sfrpermass = 1. / age massfrac0 = (sfr / mass - sfrpermass[1]) / (sfrpermass[0] - sfrpermass[1]) if finite(massfrac0) eq 0 then begin mass = mass * [massfrac0, 1-massfrac0] endif else begin massfrac1 = (sfr / mass - sfrpermass[0]) / (sfrpermass[1] - sfrpermass[0]) if finite(massfrac1) eq 0 then begin print, 'Error - cannot determine relative masses.' return endif mass = mass * [1-massfrac1, massfrac1] endelse endif ; specify sfr, mass -> calculate age. if n_elements(age) eq 0 and n_elements(mass) gt 0 then begin age = fltarr(npop) if n_elements(mass) eq 1 then mass = fltarr(npop) + mass if n_elements(texp) eq 1 then texp = fltarr(npop) + texp for p = 0, n_elements(npop)-1 do begin if texp[p] gt 0 and finite(texp) then begin age[p] = texp[p] * alog(mass / (sfr*texp) + 1) ;from: mass = sfr * texp * (exp(age/texp) - 1) endif if finite(texp[p]) eq 0 then begin age[p] = mass[p] / sfr[p] endif endfor ;if age lt 1e5 then begin ; print, 'WARNING: Age has converged to an unrealistic value, check initial conditions.' ;endif endif ; specify sfr, scalar age --> calculate scalar mass (unusual to do this if texp is specified) if n_elements(mass) eq 0 and npop eq 1 then begin if finite(texp) then mass = sfr * texp * (exp(age/texp) - 1) $ else mass = sfr * age endif ; specify vector mass, age -> no change necessary. endif else begin if n_elements(mass) eq 0 then print, 'Need to specify mass, not SFR if going the custom route' endelse ; make sure all true inputs are arrays of length npop. if n_elements(age) eq 1 then age = fltarr(npop) + age if n_elements(texp) eq 1 then texp = fltarr(npop) + texp ;if n_elements(metallicity) eq 1 then metallicity = fltarr(npop) + metallicity ;if n_elements(embedfrac) eq 1 then texp = fltarr(npop) + texp ;if n_elements(embedtime) eq 1 then embedtime = fltarr(npop) + embedtime ;if n_elements(Av) eq 1 then Av = fltarr(npop) + Av ;if n_elements(embedAv) eq 1 then embedAv = fltarr(npop) + embedAv if n_elements(mass) lt npop then begin print, 'Something has gone wrong - mass is insufficiently determined. Check inputs' stop endif ; SFR option control if (n_elements(sfr) gt 0. or n_elements(vissfr) gt 0.) and n_elements(embedfrac) eq 0 and n_elements(embedsfr) eq 0 then begin if n_elements(vissfr) eq 0 then vissfr = sfr if n_elements(sfr) eq 0 then sfr = vissfr embedfrac = 0. embedsfr = 0. endif if n_elements(sfr) gt 0 and n_elements(embedfrac) gt 0 then begin vissfr = sfr * (1-embedfrac) embedsfr = sfr * embedfrac endif if n_elements(vissfr) gt 0 and n_elements(embedfrac) gt 0 then begin sfr = vissfr / (1-embedfrac) embedsfr = sfr * embedfrac endif ; construct star formation histories if n_elements(sethist) eq 0 then begin sfhist = fltarr(221) for p = 0, npop-1 do begin if texp[p] eq 0 then sfhist += mass[p] * popsynhist(age=age[p], model='ssp', tsmooth=tsmooth) if texp[p] gt 0 and finite(texp[p]) then sfhist += mass[p] * popsynhist(age=age[p], model='exp', texp=texp[p], tsmooth=tsmooth) if finite(texp[p]) eq 0 then sfhist += mass[p] * popsynhist(age=age[p], model='cont', tsmooth=tsmooth) endfor endif else begin sfhist = sethist.sfr endelse ; in situations of truncated-to-zero recent SFR, numerical issues can cause mildly negative SFRs. ; ideally should prevent this happening further upstream, but deal with it here for now. wsomewhatnegsfr = where(sfhist lt 0 and abs(sfhist) le 2d-5*max(abs(sfhist)), ctsomewhatneg) if ctsomewhatneg gt 0 then sfhist[wsomewhatnegsfr] = 0. wnegsfr = where(sfhist lt 0, ctneg) if ctneg gt 0 then begin ;print, 'Negative SFR detected; returning null spectrum.' ; should print this if "loud" set...? Llambda = wav*0. neb = wav*0. Lyoung = wav*0. Lold = wav*0. endif else begin Llambda = popsynsed(wav, sfhist, mass=total(mass), metallicity=metallicity, Lyoung=Lyoung, Lold=Lold, $ cutage=embedtime, bct=bct, $ binning=binning, resolution=resolution) ; this is in solar luminosities/Ang, but this program works in CGS and uJy if keyword_set(nonebular) eq 0 then begin neb = calcnebular(wav, sfr=sfhist[0], resolution=resolution, gasmetal=gasmetal, q=q, $ lines=lines, lineluminosity=lineluminosity, linewav=linewav) ;linewav = [6563., 4861., 4341., 3727., 5007., 6570., 1217.] lineluminosity = redden(linewav, lineluminosity, av, rv, dustlaw=dustlaw) endif else begin neb = wav*0 endelse endelse if n_elements(lyoung) eq 0 then begin print, 'Something has gone wrong - SED synthesis failed?' return endif Lyoung += neb ; sigh, this old/young business is such a kludge. Llambda += neb hist = {t:bct, sfr:sfhist, sfrobs:sfhist*embedfrac*(bct lt embedtime)} Llambda_embed = Lyoung*embedfrac Llambda_vis = Lyoung*(1-embedfrac) + Lold nwav = n_elements(wav) ;dwav = (shift(wav,-1) - shift(wav,1))/2. ;dwav[0] = wav[1]-wav[0] ;dwav[nwav-1] = wav[nwav-1]-wav[nwav-2] dwav = diff(wav) Lbol = total(Llambda*dwav) Llambda_vis = redden(wav, Llambda_vis, av, rv, dustlaw=dustlaw) Llambda_embed = redden(wav, Llambda_embed, embedav, rv, dustlaw=dustlaw) Lvis = total(Llambda_vis*dwav) Lembedvis = total(Llambda_embed*dwav) Lir = Lbol - (Lvis+Lembedvis) ;print, 'L_bol: ', Lbol ;print, 'L_vis: ', Lvis+Lembedvis, ' (',clip(Lvis),',',clip(Lembedvis),')' ;print, 'L_IR: ', Lir ; fracabsorbed = Lir/Lbol ; not actually used anymore here ; if Lbol eq 0. then fracabsorbed = 0. ; prevent NaNs LL = (Llambda_vis+Llambda_embed) * wav ; LL = lambda*L_lambda or nu*fnu, solar luminosities ; freq = 3d+18/wav ; Hz ; should be unnecessary as this was already set Lnu_cgs = LL*3.839d+33/freq Lnu = Lnu_cgs ; freqGHz = freq / 1e+9 return end ; Primary procedure - see top of file for details. pro buildgalaxy, freq, Lnu, $ sfr=sfr, tav=tav, age=age, model=model, texp=texp, tsmooth=tsmooth, metallicity=metallicity, mass=mass, $ embedfrac=embedfrac, av=av, embedav=embedav, dustlaw=dustlaw, rv=rv, $ Tdust=Tdust, beta=beta, alpha=alpha, gasmetal=gasmetal, q=q, $ binning=binning, resolution=resolution, $ Ltot=Ltot, Lbol=Lbol, Lir=Lir, hist=hist, $ lines=lines, lineluminosity=lineluminosity, sethist=sethist, nonebular=nonebular ; note: Lbol and Lir are now ***outputs***!!! But Ltot is an input. if n_elements(Ltot) gt 0 then if Ltot gt 0 then begin if n_elements(sfr) gt 0 then if sfr gt 0 then print, 'Cannot specify both SFR and Ltot.' if n_elements(mass) gt 0 then if mass gt 0 then print, 'Cannot specify both mass and Ltot.' ; This is very weak right now. In the future, though, you should be able to do this by ; integrating the BC03 curves (ahead of time) and normalizing to luminosity within calcstellar. endif if n_elements(av) eq 0 then begin print, ' Need to specify av (dust extinction column)' return endif if n_elements(av) eq 0 then begin print, ' Need to specify mass (stellar mass)' return endif if n_elements(age) gt 2 then begin print, ' More than two components no longer supported.' return endif calcstellar, freqstell, Lnu_stellar, $ sfr=sfr, tav=tav, age=age, texp=texp, tsmooth=tsmooth, metallicity=metallicity, mass=mass, $ embedfrac=embedfrac, av=av, embedav=embedav, dustlaw=dustlaw, rv=rv, $ gasmetal=gasmetal, q=q, $ binning=binning, resolution=resolution, $ Lbol=Lbol, Lir=Lir, hist=hist, sethist=sethist, nonebular=nonebular, $ lines=lines, lineluminosity=lineluminosity if n_elements(Ltot) gt 0 then begin Lnu_stell *= Ltot/Lbol Lir = Lir * Ltot/Lbol Lbol = Ltot endif sfr0 = hist.sfr[0] ; lowest bin pahfrac = 0.1 ;if n_elements(metallicity) then if metallicity lt 0.3 then $ ; pahfrac = 0.1 * (metallicity/0.3) ; See Figure 17 of Smith et al. if n_elements(gasmetal) then if gasmetal lt 0.3 then $ pahfrac = 0.1 * (gasmetal/0.3) ; See Figure 17 of Smith et al. minstellfreq = 3e4 ; maximum frequency in GHz for using the stellar wavelength grid longfreq = (10d)^(-1 + (1+alog10(minstellfreq))*dindgen(1000)/1000.) ; covers 10^-1 to 3x10^4 GHz, or 10um to 3m ; formerly 1 GHz -> 10^16 Hz, or 0.3m -> 300 Ang (0+7*) ; in the future, should give an option to keep the output grid the same freq = [longfreq, freqstell[reverse(where(freqstell gt minstellfreq))]] wav = 3d9/freq ; angstroms presumably Lnu_stell_long = interpol(Lnu_stellar, freqstell, longfreq) Lnu_stellar = [Lnu_stell_long[where(longfreq le minstellfreq)], Lnu_stellar[reverse(where(freqstell gt minstellfreq))]] ; the optical power absorbed by dust is reradiated: 10% as PAHs, 90% as thermal continuum. ; (measured PAH fractions in reality range between 3-16%; Smith, Draine et al.) Lnu_nonthermal = calcnonthermal(freq, L=luminosity, sfr=sfr0, alpha=alpha) Lnu_thermal = calcthermal(freq, T=Tdust, L=Lir*(1-pahfrac), beta=beta) Lnu_pah = calcpah(freq, L=Lir*pahfrac) Lnu_nonstellar = Lnu_nonthermal + Lnu_thermal + Lnu_pah calz_unred2, wav, Lnu_nonstellar, -av/3. ; dust is self-absorbed at short wavelengths. Ignore dust type here. ; not sure if this is the right treatment. Note that this energy vanishes, need to feed it back in another iteration. Lnu = Lnu_stellar + Lnu_nonstellar ;Lnu = [Lnu_long[where(freq lt 10e+4)], Lnu_stellar[reverse(where(freqstell gt 10e+4))]] ; use stellar frequency indices >10^4 GHz or <30 microns, use logarithmic indices elsewhere checkplot = 0 if checkplot then begin !p.multi = 0 !p.position = 0 plot, wav/1d4, Lnu, xrange=[1e-2,1e5], yrange=max(Lnu)*[1d-12,1], /xlog, /ylog, /xsty, psym=3 oplot, 3d9/minstellfreq/1d4*[1,1], [1, 1d40], linestyle=2, color=rgb(120,120,0) oplot, wav/1d4, Lnu_nonthermal, color=rgb(255,255,0), psym=3 oplot, wav/1d4, Lnu_thermal, color=rgb(255,0,0), psym=3 oplot, wav/1d4, Lnu_pah, color=rgb(255,120,0), psym=3 oplot, wav/1d4, Lnu_stellar+2d20, color=rgb(0,0,255), psym=3 oplot, wav/1d4, Lnu, psym=3 endif ; this is no longer used ; kband = where(freq gt 3e+14/2.3/1e+9 and freq lt 3e+14/2.05/1e+9) ; slumk = total((freq[kband]-freq[kband-1]) * Lnu[kband] ) / (max(freq[kband])-min(freq[kband])) ; apparentstellarmass = 2.67d-48 * 0.4 * 1d+29 * slumk ; from Castro Ceron 2010. ; ; equivalently, slumK_sun = mag2ujy(3.28,'K')*1d-29 * (4.*3.1415 * (10*3.08d+18)^2) = 3.875d+18 ; ; stellarmass = 0.4 * slum / slumK_sun ( = 0.4 * slum * 2.58e-19, which is nearly the same) ; f = f * igm_atten(w, z) ;/(4*3.1415*(3.085d+24*dLMpc)^2) end ;Lnucgs * (1+z)^(1+alpha) / (4*3.1415*dLcm^2) * (freqGHz/1.4)^(-alpha)