[yt-users] potential unit error in cosmology.py

[Pengfei Chen]

Hi Britton,

Thanks very much for your test! The error happens when I make light rays
from one single dataset. I tried to reproduce the error with the dataset
under src/yt-hg/tests. The code I used is here:

http://paste.yt-project.org/show/6751/
After running the code, I did

h5ls -d lightray.h5/grid/redshift
and found the last redshift was 9.98999107782409. However, the correct
value should be ~ 9.979.

When I call LightRay this way, it will call the _deltaz_forward function
in cosmology_splice.py, which will call the comoving_radial_distance in
cosmology.py:

        distance2 = self.cosmology.comoving_radial_distance(z2, z)
, where distance2 is not in comoving units. Then it will calculate z2:

            z2 = ((target_distance - distance2) / m.in_units("Mpccm / h"))
+ z2

. I think the subtraction here is the problem, since it tries to convert
distance2 to comoving units, even though the numerical value of distance2
is already comoving. So a wrong factor of (1+z) is multiplied.

Any comment is appreciated.

Thanks again for your help!
Pengfei

On Mon, Aug 8, 2016 at 7:10 AM, Britton Smith <brittonsmith at gmail.com>
wrote:

> Hi Pengfei,
>
> I'm not able to reproduce what you're seeing, but before I go into that, I
> can explain a few things about how the cosmology calculator works.  In yt,
> we are only able to define a comoving unit system in relation to a proper
> unit system at a specific redshift.  We are unable to define a comoving
> unit system that is the base unit system.  This is a problem for the
> cosmology calculator, which only calculates comoving quantities.
> Therefore, in this case, we are forced to set the comoving and proper unit
> systems to be the same.  You can see this for yourself by running the
> following script:
>
> from yt.utilities.cosmology import Cosmology
> co = Cosmology()
> print (co.comoving_radial_distance(1.92, 2).to("Mpc/h"))
> print (co.comoving_radial_distance(1.92, 2).to("Mpccm/h"))
>
> Both will be the same answer.
>
> Now, back to the issue you're seeing.  I ran the following script:
> http://paste.yt-project.org/show/6741/
> to calculate the redshifts dumps needed for a 80 Mpc/h box to go from z =
> 2 to z = 1.92 and got the following result:
> CosmologyOutputRedshift[0] = 2.000
> CosmologyOutputRedshift[1] = 1.926
>
> This is pretty close to what you expected.  Additionally, if I run the
> following script,
> http://paste.yt-project.org/show/6743/
> with the above redshifts added to the parameter file, I get that dz_max
> (the change in redshift when moving 80 Mpccm/h at z = 0)
> of 0.0746521607364, which also seems to agree with what you are finding.
>
> Can you share with me the simulation parameter file and script that you're
> using to make your light rays?  I can take a look at that and see if I
> still get what I expect.
>
> Britton
>
> On Mon, Aug 8, 2016 at 6:13 AM, Pengfei Chen <madcpf at gmail.com> wrote:
>
>> Hi all,
>>
>> I was trying to generate light rays using yt 3.4-dev, but found that the
>> light rays have wrong redshift intervals. For example, my simulation box is
>> 80Mpccm/h per side, and I expect a light ray generated in a data dump at
>> z=2 with length =1 simulation unit (i.e. 80Mpccm/h) would have a redshift
>> interval from z=2 to z=1.92. However, the light ray generated by yt gives
>> spans from z=2 to z=1.97.
>>
>> As I look into this, I found that the function comoving_radial_distance
>> in cosmology.py might return the wrong unit. I think the return value
>> should be in comoving units, instead of physical units. To see this
>> directly, I made the following change in cosmology_splice.py:
>>
>> (root) ~/yt-conda/src/yt-hg/yt/analysis_modules/cosmological_observation
>> $diff cosmology_splice.py cosmology_splice.py0
>> 373d372
>> <             print target_distance, distance2, z2
>>
>> And it shows:
>>
>> 80.0 Mpccm/h 4.65320708035e+26 cm 1.97387969592 dimensionless
>> 80.0 Mpccm/h 1.19455177073e+26 cm 1.97343368065 dimensionless
>> 80.0 Mpccm/h 1.21507519065e+26 cm 1.97342592996 dimensionless
>>
>> Then I made the following change in cosmology.py:
>>
>> (root) ~/yt-conda/src/yt-hg/yt/utilities $diff cosmology.py cosmology.py0
>> 111,112c111,112
>> <         return self.quan((self.hubble_distance() *
>> <                 trapzint(self.inverse_expansion_factor, z_i,
>> z_f)).value, 'cmcm')
>> ---
>> >         return (self.hubble_distance() *
>> >                 trapzint(self.inverse_expansion_factor, z_i,
>> z_f)).in_base(self.unit_system)
>>
>>
>> Then I get:
>> 80.0 Mpccm/h 4.65320708035e+26 cmcm 1.92163908776 dimensionless
>> 80.0 Mpccm/h 3.62771515661e+26 cmcm 1.92124702027 dimensionless
>>
>> With the change of unit from cm to cmcm, the light rays have the right
>> span.
>>
>> Even though this solves my problem, I am not sure if similar problem
>> still exists. For example, instead of making change in
>> comoving_radial_distance, we might need to change hubble_distance into
>> comoving units. Hopefully someone familiar with yt unit system could check
>> this.
>>
>> Thanks,
>> Pengfei
>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>
>>
>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20160808/905d2d7f/attachment.html>