[yt-dev] bug in SFRD calculation

[Nathan Goldbaum]

Hi Geoffrey,

I'm not very familiar with the star formation rate calculator, but if you
think there is a bug I'd encourage you to submit a pull request with a fix.

-Nathan


On Thu, Apr 3, 2014 at 6:37 PM, Geoffrey So <gsiisg at gmail.com> wrote:

> There seems to be a bug in the star formation rate density calculation
> routine.  The volume is only the real 'proper volume' at the current
> redshift and not at any other redshifts, because it uses the current
> redshift in the calculation of proper volume in all redshifts.
>
> #---------------------------------
> from yt.mods import *
> from yt.analysis_modules.star_analysis.api import *
> import pylab
>
> pf=load('HD1200/HD1200')
> dd=pf.h.all_data()
> sfr = StarFormationRate(pf, data_source=dd, bins=10)
> sfr.write_out(name="StarFormationRate.out")
> time, lookback, z, Msolyr, MsolyrMpc3, Msol,
> cumMsol=na.loadtxt('StarFormationRate.out',usecols=(0,1,2,3,4,5,6),skiprows=2,unpack=True)
>
>
> cMpc3,=pylab.plot(z,Msolyr/(pf['CosmologyComovingBoxSize']/pf['CosmologyHubbleConstantNow'])**3)
> currentz,=pylab.plot(z,MsolyrMpc3/(1+pf.current_redshift)**3+1e-5)
> Expect,=pylab.plot(z,MsolyrMpc3/(1+z)**3)
>
> pylab.legend([cMpc3,currentz,Expect],
>              ['[s$^{-1}$cMpc$^{-3}$]',
>               'use only current z to fit',
>               'Expected to fit'])
>
> pylab.savefig('SFRD.eps')
> #---------------------------------
>
>
> in the file:
> yt/analysis_modules/star_analysis/sfr_spectrum.py
> line 151:
>             self.Msol_yr_vol.append(self.mass_bins[i] / \
>                 (self.time_bins_dt[i] * tc / YEAR) / vol)
>
> The vol is a fixed number instead of varying by the redshift, even though
> the numerator in solar masses per year is from different redshifts.
>
> So I think a couple lines above it should create an empty array of vol
> vol = []
> In the loop
> for i, time...
> calculate the right volume vol[i]
> Then instead of dividing by the constant vol, divide by vol[i]
> corresponding to the right redshift.
>             self.Msol_yr_vol.append(self.mass_bins[i] / \
>                 (self.time_bins_dt[i] * tc / YEAR) / vol[i])
>
> From
> G.S.
>
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-dev-spacepope.org/attachments/20140403/1d17750d/attachment.htm>